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Motivated by the gauge/gravity duality, we introduce a numerical scheme based on generalized 
harmonic evolution to solve the Einstein field equations on asymptotically anti-de Sitter (AdS) 
spacetimes. We work in global AdSs, which can be described by the (t,r,x,0,cj>) spherical coordi- 
nates adapted to the RxS 3 boundary. We focus on solutions that preserve an SO(3) symmetry that 
acts to rotate the 2-spheres parametrized by 8,(j>. In the boundary conformal field theory (CFT), 
the way in which this symmetry manifests itself hinges on the way we choose to embed Minkowski 
space in R x S 3 . We present results from an ongoing study of prompt black hole formation via 
scalar field collapse, and explore the subsequent quasi-normal ringdown. Beginning with initial 
data characterized by highly distorted apparent horizon geometries, the metrics quickly evolve, via 
quasi-normal ringdown, to equilibrium static black hole solutions at late times. The lowest angular 
number quasi-normal modes are consistent with the linear modes previously found in perturbative 
studies, whereas the higher angular modes are a combination of linear modes and of harmonics 
arising from non-linear mode-coupling. We extract the stress energy tensor of the dual CFT on 
the boundary, and find that despite being highly inhomogeneous initially, it nevertheless evolves 
from the outset in a manner that is consistent with a thermalized N = 4 SYM fluid. As a first 
step towards closer contact with relativistic heavy ion collision physics, we map this solution to a 
Minkowski piece of the RxS 3 boundary, and obtain a corresponding fluid flow in Minkowski space. 



I. INTRODUCTION 

It has been hoped since the beginnings of the 
AdS/CFT correspondence, first formulated in [IH3], that 
gauge-gravity dualities would eventually be used to re- 
late gravitational calculations to experimentally testable 
predictions in gauge theory. Over the past several 
years, there has been a flurry of activity in applying 
AdS/CFT to problems in high energy and condensed 
matter physics, which began with the discovery of the 
strongly-coupled quark-gluon plasma formed in heavy ion 
collisions. Recent proposals have suggested application 
to other systems, including those that exhibit supercon- 
ductivity, the quantum hall effect, and superfluidity. For 
some review articles see [4TU1|. Though conformal field 
theories cannot model the exact properties of all of these 
physical phenomena, it is hoped that they nevertheless 
capture aspects of the essential physics. In the strong 
coupling limit of the boundary CFT, the duality maps 
CFT states to bulk gravity configurations, and there are 
many cases where the latter theory is more tractable than 
the former. In many if not most of the model problems 
studied to date, the gravity description has involved black 
holes; aside from the interesting philosophical questions 
this poses, the implication is that one needs to study bulk 
solutions within the strong- field regime of general relativ- 
ity. In situations where exact solutions are not known, 
or perturbative expansions about known solutions are in- 
adequate to capture the non-linear dynamics, numerical 
solution of the Einstein field equations for an asymptot- 
ically AdS spacetime is required. 

Numerical relativity has seen significant progress in re- 
cent years modeling dynamical, strong-field geometries, 
though the majority of applications have been to com- 
pact object collisions in asymptotically flat spacetimes. 



For some review articles see [T^U5| . The asymptotic 
structure of AdS is drastically different from that of fiat 
space; in particular the boundary of AdS is time-like 
and in causal contact with bulk geometric structures, 
on time scales relevant to the boundary physics. This 
poses unique challenges for numerical evolution. The 
majority of existing literature on numerical solution in 
AdS has focused on black hole formation in spherically- 
symmetric 3-dimensional [UtHTB], 4-dimensional [T5], 5- 
dimensional [5U] and general D-dimensional [3TJ [52] 
spacetimes, i.e. 1 + 1 dimensional simulations where the 
boundary is a single point, significantly simplifying its 
treatment. A notable exception is a study of colliding 
Shockwaves in 5-dimensional AdS (AdSs), with applica- 
tion to heavy ion collisions [23] (and see (24] [25] for follow- 
up studies). Planar symmetry was imposed in two spa- 
tial dimensions, reducing the problem to a 2 + 1 dimen- 
sional simulation. Their approach was based on a null 
(characteristic) evolution scheme which is well-adapted 
to describing such colliding plane waves. This method 
simplifies the treatment of the AdSs boundary, though 
is difficult to generalize to situations with less symmetry. 
A more recent related study of boundary hydrodynam- 
ics via numerical solution of the full field equations in 
the bulk was presented in [55] ; though the evolution was 
effectively 1 + 1 dimensional, the space-plus-time decom- 
position, as used here, in principle allows for a straight- 
forward extension to situations with less symmetry. 

Inspired by the growing success of gauge/gravity du- 
alities, we are initiating a new program to solve the 
Einstein field equations in asymptotically AdSs space- 
times, based on the generalized harmonic (GH) evolution 
scheme presented in [27, 28 . These methods were intro- 
duced in the context of asymptotically flat spacetimes, 
and it turns out to be a rather non-trivial exercise to 
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adapt them to AdSs. The main purpose of this paper 
is to describe, in detail, the steps taken towards achiev- 
ing a stable Cauchy evolution code in spacetimes that 
are asymptotically AdSs. In this first study we impose 
an SO(3) symmetry for simplicity; since the method is 
based on Cauchy evolution, it should be straightforward 
to relax this symmetry restriction for future studies, at 
the expense of computational complexity. As a sample 
application, we study the quasi-normal ringdown of ini- 
tial data describing highly distorted black holes formed 
from scalar field collapse, and extract the stress energy 
tensor of the dual conformal field theory defined on the 
E x S 3 boundary. The 50(3) symmetry we impose here 
is precisely the one identified in analytical studies of col- 
liding shock waves in AdS 5 [5^1 [3U] . 

An outline of the rest of the paper is as follows. We 
begin in Sec. [TT] with relevant facts about asymptotically 
AdSs spacetimes. In Sec. Ill we describe a method for 



constructing AdSs initial data based on a conformal de- 
composition, restricted to time-symmetric initial condi- 
tions. For this study we use a scalar field to source devi- 
ations from pure AdSs, and find that we can solve for 
arbitrarily strong initial data: the initial slice can be 
made to contain an apparent horizon, and this horizon 
can be made arbitrarily large by adjusting (for exam- 
ple) the amplitude of the initial scalar field profile. In 
Sec. |IV| we describe the GH formalism that we use to 
evolve the initial data forward in time. Crucial to the 
stability of this scheme is the asymptotic nature of the 
so-called source functions that in GH evolution are tra- 
ditionally associated with coordinate degrees of freedom. 
Here, in contrast to asymptotically flat spacetimes, we 
find that their asymptotic form can not be freely speci- 
fied if the metric deviation is to be non-singular in the 
approach to the AdSs boundary. This and other sub- 
tleties related to evolution are clarified and addressed in 
the latter parts of the section. We discretize the field 
equations using finite difference methods, and the details 
of the numerical methods that we use to solve the initial 
data and evolution equations are discussed in Sec. [V] 

Results from a study of highly deformed black holes, 
their subsequent evolution and ringdown, and the stress 
tensor of the corresponding states in the dual boundary 
CFT, are presented in Sec. |VI| What is novel about this 
study is that the initial horizon geometry cannot be con- 
sidered a small perturbation of the final static horizon, 
and hence we are probing an initial non-linear phase of 
the evolution of the bulk spacetime. Shortly after the 
initial time, the metric can be described by a combina- 
tion of quasi-normal modes and what appear to be gauge 
modes. We find frequencies that are consistent with the 
linear modes found in perturbative studies of static black 
holes, and in modes at higher angular number, we find 
evidence of non-linear mode-coupling. On the bound- 
ary, the dual CFT stress tensor behaves like that of a 
thermalized Af = 4 super- Yang-Mills (SYM) fluid. We 
find that the equation of state e = 3P of this fluid is 
consistent with conformal invariance (here, e and P are 



the rest frame density and pressure of the fluid, respec- 
tively), and that its transport coefficients match those 
previously calculated for an Af = 4 SYM fluid via holo- 
graphic methods. Modulo a brief transient that is numer- 
ical in nature, this matching appears to hold from t = 
onwards. This is not a priori inconsistent with earlier 
results from quasi-normal modes |31H33j and numerical 
analyses 23-26 , where a certain "thermalization time" 
is observed before the boundary physics begins to cor- 
respond to that of a fluid. In contrast to those studies, 
we begin with a large, distorted black hole, i.e. an inho- 
mogeneous thermal state, and the physics studied here is 
more of equilibration rather than of thermalization. 

In the final section, we transform solutions computed 
in global AdS onto a Minkowski piece of the bound- 
ary, and examine the temperature of the corresponding 
fluid flows. Under this transformation, the spatial profile 
of temperature at the initial time resembles a Lorcntz- 
flattened pancake centered at the origin of Minkowski 
space. By interpreting the direction along which the data 
is flattened as the beam-line direction, our initial data 
can be thought of as approximating a head-on heavy ion 
collision at its moment of impact. 

We relegate to Appendix |A"| some technical details on 
the effect of matter backreaction on the asymptotic met- 
ric fall-off. In Appendix [Bj we discuss how our asymp- 
totic metric boundary conditions relate to the bound- 
ary conditions derived in [31] for linearized gravitational 
perturbations on an AdS background that lead to well- 
defined dynamics, and in Appendix [C] we tabulate the 
scalar field quasi-normal mode frequencies. Throughout, 
we use geometric units where Newton's constant G and 
the speed of light c are set to 1. 



II. GRAVITY IN ASYMPTOTICALLY ADS 5 
SPACETIMES 

A. The AdS 5 Spacetime 

The Lagrangian density of gravity with cosmological 
constant A, coupled to matter with a Lagrangian density 
C m , is given by 



C = (R - 2A) + Cr. 

107T 



(1) 



where R is the Ricci scalar. The corresponding equations 
of motion then take the local form 



8ttT, 
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(2) 



Here, R^ v is the Ricci tensor, is the metric tensor and 
T^ v is the stress energy tensor of the matter. The met- 
ric of AdSs is the maximally symmetric vacuum solution 
of ([2]) in D = 5 dimensions, with a negative cosmological 
constant A = A 5 < 0. In terms of the global coordinates 
x 1 - 1 = (t,r,x,0,4>) that cover the whole spacetime, this 
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solution is 



functions 



ds 2 



g^dx^dx" 
-f(r)dt 2 + 



(3) 



./'(< 



-dr 2 



2 for convenience. 
2 ) is the metric 



where we have defined f(r) = 1 + r /L 
Here, dQ 3 2 = d\ 2 + sin 2 x {dO 2 + sin 2 
of the 3-sphere parametrized by angles x, 9, </>, and L is 
the AdS radius of curvature, related to the cosmological 
constant by A 5 = — (D - 1)(D - 2)/{2L 2 ) = -6/L 2 . We 
are free to set the AdS length scale L, which is the scale 
with respect to which all other lengths are measured. In 
the code we set L = 1, though we will continue to ex- 
plicitly display L in all of the following, unless otherwise 
indicated. 

Due to the significance of the AdSs boundary in the 
context of AdS/CFT, it is useful to introduce a "com- 
pactified" radial coordinate p so that the boundary is at 
a finite p. We choose 



(4) 



1-p/V 



where i is an arbitrary compactification scale, indepen- 
dent of the AdS length scale L, such that the AdSs 
boundary is reached when p = I. In our code and in all 
of the following, we set I = 1, though note that this scale 
is implicitly present since p has dimensions of length. 
Transforming to the p coordinate, the metric ^ takes 
the form: 



ds 2 = 



1 



(1-P) S 



-f(p)dt 2 



1 



f(p) 



dp 2 + p 2 dfl 3 



(5) 



where f{p) = (1 - p) 2 + p 2 /L 2 . 

AdS 5 can also be described as the universal cover of 
the hyperboloid X A in K 4,2 defined by the locus 



i=4 



(A- 1 ) 2 -(A°) 2 + ^(A 



i\2 



(6) 
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This space has symmetry group 50(4, 2) whose trans- 
formations preserve the quadratic form (j6j). The metric 
of AdSs is then the metric g^ v induced on the above 
hyperboloid from the flat metric Gab of an K 4,2 am- 
bient space. Here the metric Gab is simply given by 
G = diag(— 1, — 1, 1, 1, 1, 1). The hyperboloid can be 
more efficiently described in terms of embedding coordi- 
nates x^ defined by a set of embedding functions X A (x 11 ), 
so that the induced metric <7 M „ on the hyperboloid is 



fdX J 



V dx^ 



dX 1 

dx v 



G 



AB- 



(7) 



The global coordinates x M = (t, r, x, 9, 4>) can then be 
thought of as corresponding to a choice of embedding 



x- 1 


= Vr 2 + L 2 


cos(t/L) 


x° 


= \/r 2 + L 2 


sin(i/L) 


X 1 


= r sin x sin i 
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X 2 


= r sin x sin i 


9cos</> 


X 3 


= r sin x cos 
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X 4 


= r cos x- 





(8) 



Notice that t = and t = 2irL are identified on the 
hyperboloid, so that there are closed time-like curves in 
this space. AdSs is defined as the hyperboloid's universal 
cover precisely to remove these closed time-like curves. 
This universal cover is obtained by unwrapping the S 1 
parametrized by t on the hyperboloid, which would then 
run from — oo to oo in AdSs. 



B. Boundary of AdS 5 

The boundary of AdSs differs dramatically from that 
of asymptotically flat spacetimes. To see this, let us first 
remember how the analogous story unfolds in the familiar 
setting of Minkowski space, whose metric in polar coordi- 
nates ij l j iU dx fi dx L ' = —dt 2 + dr 2 + r 2 dVl 2 can be rewritten 
through a series of coordinate transformations u± — t±r, 
u ± = arctanu-t = (T ± R)/2 to read 



ds 2 



1 



4 cos 2 Uj- cos 2 u_ 



(-dT 2 + dR 2 + sin 2 Rdn 2 ) 



(9) 

Through this conformal compactification, the infinite re- 
gion t G (— oo, oo), r € [0, oo) is mapped to the interior of 
a compact region defined by \T ± R\ < it. Given that T 
is a time-like coordinate and R > is a space-like coor- 
dinate, this region is the triangle T £ [—71", it], R £ [0, w] 
in the T, R plane — see Fig. [l] Consequently, the bound- 
ary consists of the two null surfaces of future/past null 
infinity, meeting at spatial infinity. Spatial infinity in 
Minkowski space is thus not in causal contact with its 
interior. Notice that the geometry of Minkowski space is 
conformal to a patch T £ [—%, tt] of the Einstein static 
universe. 

The conformal compactification of AdSs is achieved 1 
by the spatial coordinate transformation r/L = tani?, 
which brings its metric ^ into the form 



ds 2 



cos 2 R 



^dt 2 



dR 2 + sin 2 Rdn 3 2 ) 



(10) 



where the infinite region r £ [0, oo] is mapped to the 
interior of a compact region R £ [0,n/2] — see Fig. [2] 



1 In practice, we compactify using Q, but for the current didac- 
tic discussion we choose the compactification that most closely 
resembles its Minkowski space analog. 
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FIG. 1: The conformal diagram of Minkowski space. The 
boundary consists of the point at spatial infinity i°, and 
the null surfaces at future null infinity J + and past null 
infinity J~ . In this compactification, future time-like in- 
finity i + and past time-like infinity i~ are represented by 
points. Dashed lines are constant t surfaces, and solid lines 
are constant r surfaces. 



FIG. 2: The conformal diagram of anti-de Sitter space. The 
boundary consists of the time-like surface J\ past and fu- 
ture time-like infinity are represented by the points i~ and 
i + , respectively. Conventions used here are the same as 
those in Fig. [T] 



C. Asymptotically AdSs Spacetimes 



This halved range of R implies that anti-de Sitter space 
is conformal to one-half of the Einstein static universe. 
More crucially, we have not rescaled the time coordinate 
t in this compactification. The important consequence 
is a spatial infinity that runs along the t direction: it is 
time-like and thus causally connected to the interior. 

A proper treatment of the AdS boundary is crucial to 
a solution of the Cauchy problem in an asymptotically 
AdS spacetime. Without some specification of bound- 
ary conditions at time-like infinity, only a small wedge 
to the causal future of an initial space-like foliation can 
be solved for. This is in contrast with asymptotically 
Minkowski spacetimes, where the specification of initial 
data on a slice that reaches spatial infinity is sufficient 
to evolve the entire interior. Such an approach to initial 
data is not useful in the asymptotically AdS case, par- 
ticularly in problems that are relevant to AdS/CFT: for 
these problems, the time-like boundary must be included 
as part of the spacetime. 



An asymptotically AdSs (AAdSs) spacetime is one 
that shares the boundary structure of AdSs. In partic- 
ular, it must have 50(4,2) as its asymptotic symmetry 
group. To write down the asymptotics of the fields in 
such a spacetime, let us decompose the metric by writing 



Q^iv "F h^iv 



(11) 



where is the deviation of the full metric from 
the metric of AdSs ■ For an explicit form of the met- 
ric that we evolve, the reader may wish to briefly skip 
to ((54 1. 



The matter-free asymptotics of h^ v corresponding to 
an AAdSs spacetime were found in [35], by requiring 
that these deviations h^ v satisfy boundary conditions 
at spatial infinity that are invariant under the SO (A, 2) 
symmetry group of AdSs. The resulting boundary con- 
ditions are consistent with the Fcffcrman-Graham con- 
struction [35] (see Appendix (b| . To obtain these bound- 
ary conditions, ;'i. r >] begins by noting that a typical gen- 
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erator £ M of 50(4, 2) has asymptotics 



takes the form 



r 
r 

c" 1 
s , r 

S.r 



0(1) 

0(r) 
0(1), 



(12) 



where the index to denotes the non-radial coordinates 
(t,x,6,<t>)- Since we need to ensure that 50(4,2) is an 
asymptotic symmetry, we are interested in satisfying the 
Killing equation in an asymptotic sense 



(13) 



so that for all generators £ M of 50(4, 2), the Lie derivative 
of each metric component along £ M approaches zero with 
the appropriate power of r as r oo . In other words, we 
are looking for an asymptotic form of the metric devia- 
tion h^ v ~ r p ^" , such that this asymptotic form is pre- 
served by the coordinate transformations that correspond 
to the generators £ M . Direct calculation reveals that (13) 
holds when p rr = —D — 1, p rm = —D, p mn = —D + 3. 
Since we are interested in the case of D = 5 dimensions, 
the vacuum boundary conditions read 

Kr = frr(t,X,6,4>)^+0(r- 7 ) 



h r 

tin 



1 



1 



+ 0(r- b ) 



/ m „(t,X,^,0)^+O(r- J ) 



(14) 



These boundary conditions are valid for vacuum 
AAdSs spacetimes, and as we shall see in the next sec- 
tion, continue to hold for spacetimes containing localized 
matter distributions with sufficiently rapid fall-off near 
the boundary. For more general matter distributions, 
fall-off conditions can be obtained via different methods; 
see for example [37J EE] • The most general set of fall-off 
conditions were found in [39] : these continue to hold for 
spacetimes with fewer restrictions on boundary confor- 
mal structure and boundary topology. 



D. Scalar Fields in Asymptotically AdSg 
Spacetimes 



Or 



Or 



to 



(16) 



rD D-2 , p 
L l r 

The static ansatz <f>(r) ~ r -A yields a quadratic equation 
for A whose solutions are the powers of allowed fall-offs. 
The asymptotics for a scalar field of mass m are thus 
described by 



.4 



B 



A = 




(17) 



(18) 



i 2 in D = d + 1 = 5 di- 



which gives A = 2 
mensions. The usual vacuum/localized-matter boundary 
conditions given in ( 14 1 can accomodate a scalar field 
with a vanishing A branch, ^4 = 0. We show in Appendix 
|A| th at this case does not alter the boundary conditions 
(14), because the scalar field falls off sufficiently quickly 



near the boundary. This also suggests why gravitational 
wave perturbations are consistent with ( 14 ) , since their 
propagation characteristics are "morally" equivalent to 
that of a massless scalar field, i.e. having A = 4 for B 
branch solutions. 

The situation becomes richer for nonlocalized-matter 
boundary conditions with a non-zero A branch, A =/= 0. 
In this case, the backreaction of the scalar field onto 
the geometry begins to alter the metric asymptotics. 
Scenarios studied in this paper require only the usual 



vacuum/localized-matter boundary conditions (14); sys- 



tems which require nonlocalized-matter boundary condi- 
tions will be addressed in a future work. 



E. The Stress Energy Tensor of the Dual CFT 

Here we briefly mention the most relevant entry in the 
AdS / CFT dictionary for this study, namely the one that 
enables us to extract the expectation value (T MJ/ ) CFT of 
the CFT stress energy tensor from the asymptotic be- 
havior of the metric: 2 



IT ) — lim —^T 



(19) 



As we will be coupling matter to gravity, it is impor- 
tant to know how the presence of matter alters the vac- 
uum boundary conditions (14). To understand just the 
asymptotics, it suffices to consider a static spherically 
symmetric scalar <fr — (f>(r) of mass to, for which the 
Klein-Gordon equation 



□ = TO 



(15) 



2 The full expression for the boundary stress tensor includes a fac- 
tor of 1/G, where G is Newton's constant, corresponding to a 
scaling as TV 2 in the large N gauge theory dual to AdSs. We 
have omitted the factor of 1/G from | |19| in keeping with our 
standard convention that G = 1. When quoting explicit numeri- 
cal results for the boundary stress tensor, we will also set L = 1. 
By doing so we are restricting to a specific count of degrees of 
freedom in the boundary theory; but since the scaling of (T^ v ) 
is straightforward, there is no significant loss of generality. 
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Here, q = 1 — p (though more generally q is a smooth 
positive scalar with a simple zero at the AAdS bound- 
ary), and ^T M „ is the Brown- York quasi-local stress ten- 
sor [ID]. For AAdSs spacetimes, the quasi-local stress 
tensor defined on a q = const, time- like hypersurface dM q 
was constructed in (41] . and is given by 



(q) T o = J_ 



L 



Here, ('>9,. 



— -S" M E /3 1/ V( Q 5' / g) is 



(20) 

the extrinsic cur- 
vature of the time-like surface dM q , S* M is a space-like, 
outward pointing unit vector normal to the surface dM q , 
E^ = — Sfj,Su is the induced 4-metric on dM q , V Q is 
the covariant derivative operator, and ^Gn„ is the Ein- 
stein tensor associated with S^ w . The last two terms in 



( 20 1 are counterterms designed to exactly cancel the di- 



vergent boundary behavior of the first two terms of ( 20 ) 
evaluated in pure AdSs. 

A feature of the stress tensor (20) is that it is non- 



vanishing even when the geometry is that of pure AdSs. 
This non- vanishing piece was already noticed in [4"T] . 
and was correctly identified as the contribution from the 
Casimir energy of the boundary CFT: this CFT is de- 
fined on a manifold with topology R x S 3 , and so can 
have a non-vanishing vacuum energy. Since this Casimir 
contribution is non-dynamical, we consider it as part of 



our background vacuum and simply subtract it from ( 20 ), 
obtaining 



{i)t — (<j)t° 



(?) 



ft V • 



(21) 



Setting L = 1, the non-zero components of the 
Casimir contribution are = 3g 2 /(647r), 

(?)+ 

{q) ta, SOT 



<?7(64^), Mtge = g 2 sin 2 x/(647r), and Wt w = 



III. INITIAL DATA 

Initial data for Cauchy evolution of the Einstein field 
equations ^ is not freely specifiable, but is subject to 
a set of D constraint equations: the Hamiltonian con- 
straint and the D — 1 non-trivial components of the mo- 
mentum constraints. There are many conceivable ways 
of finding initial data that are consistent with these con- 
straints; here, we adapt to AAdSs spacetimes the tradi- 
tional Arnowitt-Deser-Misner (ADM) [32]-based confor- 
mal decomposition approach often used in asymptotically 
flat spacetimes. See |43j for a recent review. To simplify 
this first study, we restrict initial data to a moment of 
time symmetry, and we use a scalar field to source non- 
trivial deviations from pure vacuum AdSs. 

Specifying time symmetric initial data is equivalent to 
demanding that the extrinsic curvature of the initial t = 
slice Sj=o vanishes. The extrinsic curvature of a constant 



t slice Tit of the spacetime is defined as 



where 



2 ^nlfiv j 



-ad^t 



(22) 



(23) 



is the unit time-like one-form normal to the E 4 slice, a 
is the so-called lapse function, and 7^ is the 4-metric 
induced on E t , expressible in local coordinates as 

l^v = g^v + n^riy. (24) 

The momentum constraints are the D — 1 equations 

D V K^ - ^ U D V K = 8tt^, (25) 

where = 7 A1 l/ V !y is the derivative operator intrinsic to 
E t , and 



3 



-T Q pn a r 



(26) 



is the momentum of any matter in the spacetime. Time 
symmetry K^\ t= o — requires that in order for (25) to 



be satisfied, the momentum density j^ 1 must vanish every- 
where on S t= o- This leaves us with only the Hamiltonian 
constraint; the solution of this equation is the subject of 
the next few sections. 



A. Hamiltonian Constraint 

The Hamiltonian constraint is a single equation 

(4) i? + K 2 - K^K" V - 2A 5 = MirpE, (27) 

where is the Rlcci scalar of the geometry intrinsic 
to the slice E t) and 



Pe = T iiU n tJ, n v 



(28) 



is the energy density on the slice. At a moment of time 
symmetry ( 27 1 simplifies to 



(4) i?- 2A 5 = 16tt Pe . 



(29) 



Following the conformal approach, we will solve this 
equation by requiring that our spatial 4-metric 7^ be 
conformal to the 4-metric 7 M „ of a spatial slice of vac- 
uum AdSs 

= £"2y^ (3 ) 

for some positive smooth function £ on E t=0 with bound- 



ary condition CI as = The Hamiltonian constraint (29) 
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can then be expressed as 

^R - eCH^AAC - C 2 2A 5 - 16ttC 2 P e, (31) 

where D a is the covariant derivative operator compatible 
with 7^, and ( 4 >R is the corresponding Ricci scalar. 
is readilv computed from the spatial part of the AdSs 
metric (jS}, giving a constant ^R = -12/i 2 = 2A 5 . 
This lets us rewrite (31 ) as 



1 



1 



T ? D a D Q - ^A 5 C + 5 (A 5 + 87rp B ) C a 



(32) 



Notice that the Hamiltonian constraint ( 32 ) does not con- 



tain the lapse function a; this is consistent with the un- 
derstanding that unlike a specification of the lapse a, 
which encodes the manner in which data evolves away 
from the initial slice, the Hamiltonian constraint may 
only set data intrinsic to the initial slice itself. Further- 
more, by restricting our attention to conformally AdS 
initial data, we have written this Hamiltonian constraint 
as a non-linear elliptic equation for the conformal factor 
(. This equation can be solved once we specify the mat- 
ter energy density pe, which we discuss in the following 
section. 



B. Hamiltonian Constraint with Scalar Matter 



The Hamiltonian constraint thus takes the form 



7 W £*AC - g(A 5 - torfldiWrfX 



■- (A 5 + 8ttV(£)K 3 



0. 



(36) 



The choice of <fi on the spatial slice is completely arbi- 
trary. For the tests and quasi-normal mode study de- 
scribed here, we restrict ourselves to free, massless fields 
i.e. V(<fr) = 0. For the spatial profile of these fields, 
we use the following 5-parameter generalized Gaussian 
function: 



<l>(p,x) = <7 4 (1 + p) 4 A,exp 



(R(p, X )-Ro) 
S 2 



(37) 



where 



R(p,x) 
y(p,x) 



lx(p,x) 2 , y(p,x) 2 



pcos\ 
psinx. 



(38) 



Here, Aq is the maximum amplitude, Rq fixes the radial 
position of the maximum, <5 sets the overall compactness 
of the profile, and w x ,w v can be adjusted to set the rela- 
tive compactness of the profile in the x,y directions. The 
q 4 factor ensures that this profile has the correct fall-off 
for a massless scalar field, consistent with ( 17 ) and ( 18 ) 3 ; 
this is supplemented by a (1 + p) 4 factor to maintain the 
original Gaussian profile's even character near the origin. 



In this section, we write down an explicit form for the 



matter source term in ( 32 ) , constructed from scalar field 



initial data. Scalar fields are particularly convenient for 
our purposes, since their energy density constitutes a pa- 
rameter with which we can tune the initial data. Con- 
sidering cases of incrementally larger energy density then 
allows us to approach the dynamical, strong field regime 
in a controlled fashion. 

The Lagrangian density of a scalar field <p with a po- 
tential V{<f)) is given by 

Ct4>,d tl <i>) = ~g a f } d a <i>d fs cj>-V{cl>). (33) 



Varying the action constructed from ( 33 ) with respect to 
the metric g a p gives an energy-momentum tensor 

Tftu = d^d v <\> - {^g af) d a ^d^ + V(<j>)\ . (34) 



Substituting (23), (24), (30 1, and (34 1 into (28), then us- 



ing the restriction of time-symmetry, which for the scalar 
field amounts to setting dt4>\t=a — 0, we obtain 



Choice of Coordinates in Terms of Lapse and 
Shift on the Initial Slice 



To complete the specification of initial data, we must 
choose coordinates on the initial slice, i.e. the form in 
which we wish to represent our metric 7^. This spec- 
ification fixes the remaining coordinate degrees of free- 
dom, and within the ADM decomposition, it amounts 
to a choice of the lapse function a and shift vector 
defined via 

g^dx"dx u = -a 2 dt 2 +j ij (dx i + P'dt^dx 1 + & dt). (39) 

Here, the latin indices i,j only sum over the spatial co- 
ordinates. 

For spatial coordinates, we choose the compactified 
x% — {pi X-, $1 4>) as defined in the discussion preceding dSJ). 
To understand our choice of initial shift vector f3 % and 
an initial lapse function a, we must first understand a 



PE = C 2 f j d i <Pd j ct> + V(<t>). 



(35) 



3 The q 4 prefactor on the right side of j37| l means that we are 
not deforming the CFT by the dimension 4 operator dual to the 
massless scalar d>. 



few subtleties concerning gauge choices in asymptotically 
AdS spacetimes, and so we defer this discussion until 
Sec. |IVD| In this later section, we will recast our choice 
of /3 l ,a in terms of an equivalent choice of gt v , dtgtv at 
t = 0. 



H": 



IV. EVOLUTION 



In this section we describe the ingredients we have 
found necessary to achieve stable evolution of AAdSs 
spacetimes within the GH formalism. The solutions in 
this initial study preserve an SO(3) symmetry, which is 
sufficiently general to be physically relevant, as well as 
capture many of the problems and issues that need to 
be resolved for stable evolution. Chief among these are 
(a) decomposing the metric into a form that analytically 
factors out the AdS divergences, and dividing out suffi- 
cient powers of 1 — p from what remains, allowing us to 
set boundary conditions that impose the desired leading- 
order deviation from AdS, and (b) imposing an asymp- 
totic gauge in terms of the source functions H" of the 
GH formalism that is consistent with the desired fall-off. 
These two issues are in fact intimately related in AAdS 
spacetimes — asking for coordinates where the leading- 
order metric deviations are non-singular in the approach 
to the boundary, together with a choice of the form of 
the background (singular) AdS metric, completely fixes 
any residual gauge freedom. 

The next few sections are structured as follows: in 
Sec. | IV A| we review the GH formalism, in Sec. |IVB| we 
describe the metric decomposition that we use, and in 
Sec. |IV C| we look at the asymptotic form of the field 
equations and derive the relevant gauge conditions for 
GH evolution. In Sec. |IVD| wc show how to choose 
3tv\t=oi ®t§tv\t=o 011 ^ ne initial time slice, such that they 
are compatible with these gauge conditions. 



Ox" = 



i—g dx 



- (V=99 a,t ) = -9 aP Kp = (40) 



where g is the determinant of the metric, and T^a are the 
Christoffel symbols. To see why this has proven to be so 
useful for Cauchy evolution, let us begin by rewriting the 
field equations ^ in trace-reversed form 



11V 7 



where 



+ 8tt I T^ u — -T a a g^ v 



(41) 



(42) 



When viewed as a set of second-order differential equa- 
tions for the metric g^ u , the field equations in the form 



(41| do not have any well-defined mathematical charac- 



ter (namely hyperbolic, elliptic or parabolic), and in fact 
are ill-posed. Fixing this character requires choosing a 
coordinate system. A well-known way to arrive at a set 
of strongly hyperbolic equations is to impose harmonic 
coordinates, namely (40) with H" = 0. Specifically, this 



condition (and its gradient) can be substituted into the 
field equations to yield a wave equation for the principal 
part of each metric element, g a ^d a d/3g flL , + ... = 0, where 
the ellipses denote lower order terms. 

One potential problem with harmonic coordinates, in 
particular in a highly dynamical, strong-field spacetime 
evolved via a Cauchy scheme, is that beginning from a 
well-defined initial data surface t — const, which is every- 
where space-like, there is no guarantee that t, subject to 
the harmonic condition, will remain time-like throughout 
the spacetime as evolution proceeds. If t becomes null or 
space- like at a point, standard numerical techniques will 
break down. A solution to this, first suggested in [UJ, is 
to make use of source functions (as originally introduced 
in [46; ) . Note that any spacetime in any coordinate sys- 
tem can be written in GH form; the corresponding source 
functions arc simply obtained by evaluating the defini- 
tion |40|. Thus, trivially, if there is a well-behaved, non- 
singular coordinate chart that covers a given spacetime, 
then there is a GH description of it. The difficulty in a 
Cauchy evolution is that this chart is not known a-priori, 
and the source functions H" must be treated as indepen- 
dent dynamical fields. Finding a well-behaved coordinate 



A. The Generalized Harmonic Formalism 



Here we give an overview of the generalized harmonic 
(GH) formalism with constraint damping; for more de- 
tails see [El 123 HH SI]. The GH formalism is based 
on coordinates x" that are chosen so that each coordi- 
nate satisfies a scalar wave equation with source function 



4 As can be seen from j40| l H 11 is not a vector in the sense of its 
properties under a coordinate transformation, rather it trans- 
forms as the trace of the metric connection. One can introduce 
additional geometric structure in the form of a background met- 
ric and connection to write the GH formalism in terms of "stan- 
dard" tensorial objects. However, in a numerical evolution one 
must always choose a concrete coordinate system, and hence the 
resulting equations that are eventually discretized are the same 
regardless of the extra mathematical structure introduced at the 
formal level. 
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chart then amounts to supplementing the Einstein field 
equations with a set of evolution equations for iJ M , which 
can now be considered to encode the coordinate degrees 
of freedom in our description of the spacetime. 

The field equations in GH form thus consist of the 



Einstein equations (41 1, brought into hyperbolic form via 
the imposition of ( 40 ) 



1 aB aB 



(43) 



together with the relevant evolution equations for the 
matter (here the Klein-Gordon equation (15)) and a set 
of equations for the source functions, which wc write sym- 
bolically as 



(44) 



Even though are now treated as independent func- 
tions, we are only interested in the subset of solutions 



constraints (40). Introducing 



to the expanded system (43 1,(44) that satisify the GH 



(45) 



we thus seek solutions to (43 1,(44) for which C M = 0. An 



equivalent way of obtaining (|43|) from (|41|) is to subtract 
V(^C„) from so that 



T,, 



Q. 



(46) 



The effect of this subtraction becomes obvious when 
we rewrite the Ricci tensor explicitly in terms of Ox^ 
so R n 



l„af3 

2" 9ia>,a8 



„aB 

9 ,(p9i>)a,3 



T a p^Y^ au . We see that the subtraction of V( M C„) is 
simply designed to replace the V^Da:,,) term in R^ v by 
an equivalent V^ifj,) term. We also see that a solution 
of the Einstein field equations (41) is also a solution of 
(46), as long as the constraints C M = are satisfied. 



For a Cauchy evolution of the system (43), (44), we 
need to specify initial data in the form 



with (48). Furthermore, using a contraction of the sec- 
ond Bianchi identity V M i? A1 ^ = W v R/2, one can show 
that satisfies the following hyperbolic equation: 



/ 1 V • 



(49) 



Thus, if we imagine (analytically) solving (43), (44) using 



initial data satisfying H8h supplemented wit 



Doundary 



conditions consistent with C 7 ' = on the boundary for 



all time, then ( 49 1 implies that C M will remain zero in 



the interior for all time. 



At the level of the discretized equations, however, 
is only zero up to truncation error. This is not a priori 
problematic: numerically one only ever gets a solution 
approximating the continuum solution to within trunca- 
tion error. However, experience with asymptotically-flat 
simulations suggest that in some strong-field spacetimes, 
equation ( 49 1 for C M admits exponentially growing so- 
lutions (the so-called "constraint-violating modes"). At 



any practical resolution, this severely limits the amount 
of physical time for which an accurate solution to the de- 
sired C 7 ' = Einstein equations can be obtained. In 
asymptotically flat spacetimes, supplementing the GH 
harmonic equations with constraint- damping terms as 
introduced in [47] suppresses these unwanted solutions. 
Anticipating similar problems in AAdS spacetimes, and 
that constraint damping will similarly help, we add the 
same terms to ( 43 ) , and arrive at the final form of our 



evolution equations 



1 



aB aB 
-j9 9nv,aP - 9 ,{ M 9v)a,P 

^-{il.f) ^flT [LIS r PfjY^ ' OLU 

T^Kg^v + 8tt ( T^ v - ^T a a g^ 



(50) 



Here, the unit time- like one- form is defined as in ( 23 1 , 
and the constraint damping parameters k € (—00, 0] and 
P G [—1,0] are arbitrary constants. In all simulations 
described here, we use K = — 10 and P = — 1. 5 



<?^|t=0j dt9vw\t=Oi 
subject to the constraints 

C%= = 0, 3 4 C| t=0 = 0. 



(47) 



(48) 



One can show (see for e.g. [IS]) that if (48) is satis- 
fied, then the ADM Hamiltonian and momentum con- 
straints will be satisfied at t = 0. Conversely, if the 
ADM constraints are satisfied at t = together with 
C^\t = a = (this latter condition is satisfied trivially 
computing fP| t=0 by substituting (47) into (40)), then 
9 t C M | t= o = 0. Thus, our initial data method described 
in the previous section will produce data consistent 



Note that the new terms are homogeneous in C^, and 
hence do not alter any of the properties discussed above 
that relate solutions of the Einstein evolution equations 
and ADM constraints with those of the corresponding 
GH equations, with the exception that the constraint 
propagation equation (49) picks up additional terms, 
again homogeneous in C M (see for example |4"T]). 



5 We did not perform any systematic survey by varying k or P, 
though with a little experimentation found the exact value of K 
was not too important to achieve effective constraint damping, 
though wc found that it was important to keep P close to —1. 
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B. Evolution Variables and Boundary Conditions 

The boundary is crucial for evolution in asymptotically 
AdSs spacetimes. To find the most natural variables to 
evolve in this setting, we first need to gain some intuition 
on how the fields behave near the boundary at p = 1. To 
begin, let us again use (11 1 to decompose the metric g^ v 
into a pure AdSs piece g^ and a deviation h^. From 
the point of view of the evolution equations ( 50 ) , this 
decomposition allows us to analytically eliminate a subset 
of asymptotically singular terms representing the pure 
AdSs solution. From the point of view of the boundary 



conditions (14 1, this decomposition guides our choice of 
variables that are most suitable for Cauchy evolution. 

Our choice of evolution variables is largely motivated 
by considerations similar to those that arise when nu- 
merically solving hyperbolic partial differential equa- 
tions (PDEs) whose domain includes a formally singular 
boundary, where one seeks solutions that remain regular 
at the boundary (see for example [35] , where this method 
was introduced to study the evolution of gravitational 
waves in axisymmetry in a domain that includes the axis 
of symmetry) . Before listing the variables that we use to 
represent the metric deviation in AAdS spacetime, let us 
first colloquially describe the reasoning behind their def- 
initions. We will not prove that the following is a correct 
(or complete) characterization of the AAdS boundary 



behavior of our coupled system of PDEs (15), (44), (50 1; 
rather we will take the empirical approach that if by using 
this regularization scheme we are able to obtain stable, 
convergent numerical solutions, then this strongly sug- 
gests that the regularization is consistent, at least for the 
set of initial data considered. Though note that in [34] 
a complete characterization of boundary conditions con- 
sistent with stable evolution for linearized gravitational 
perturbations on an AdS background was given. As dis- 
cussed in Appendix [B] the boundary conditions we de- 
scribe below are consistent with the Friedrich self-adjoint 
extension of the operator describing the scalar sector of 
gravitational perturbations. Thus, insofar as the linear 
problem guides the full non-linear problem, we can have 
some confidence that the following prescription is well- 
posed. 

To illustrate, consider a function fit, q, x, 9, 0) that can 
be expanded in a power series in q, where the boundary 
is located at q = 0: 



/ = fa + fiq + h<f + hq 3 + 



(51) 



Here, /o, /1, /a, /3, ■•■ are functions of t, x, 9, <j>. Now sup- 
pose that for a regular solution (which in our case means 
a solution consistent with the desired AAdS fall-off) the 
first n terms of the RHS are required to be zero, and 
that the n-plus-first term describes the leading-order be- 
havior for the particular physical solution of interest. At 
first glance, this would suggest that we need to supply 
n + 1 boundary conditions. However, if / satisfies a hy- 
perbolic PDE with "standard" characteristic structure 



i.e. with an inward and outward propagating mode in 
each spatial direction, then we are usually only free to 
choose 1 boundary condition, effectively fixing the mode 
propagating into the domain. Furthermore, at a singu- 
lar boundary where we demand regularity, one is often 
not free to choose even this ingoing mode — the outgo- 
ing mode together with regularity completely fixes it. As 
proposed in [JSJ, a solution is to define a new evolution 
variable / via 



f{t,q, X ,9,<jy) = f(t, q ,x,9A)cr 



(52) 



and demand that / satisfy a Dirichlet condition f(t, q = 
0, x, 9, 4>) — at the boundary. Plugging this into (51) 
gives 



/ = ••• + /n-2<? 1 + fn-1 + fnq + fn+iq 2 + 



(53) 



One can see that if we choose regular initial data for / 
that has f(t = 0, q — 0) = 0, this eliminates all the com- 
ponents of / with undesired fall-off at t = 0. This will 
give a regular solution (analytically) for all time if the 
differential equation for / admits a unique solution con- 
sistent with the desired fall-off. Note that by factoring 
out n — 1 powers of q, we have assumed that the nature 
of the boundary is such that we are not free to spec- 
ify the leading-order behavior encoded in f n as a formal 
boundary condition, i.e. initial data together with evo- 
lution should uniquely determine it (if this were not the 
case, one could factor out an additional power of q, and 
explicity set f n ). Again, as described in more detail in 
Appendix |B this is consistent with the analysis of [M] 
on admissib e boundary conditions for linearized grav- 
itational perturbations of AdS. This is also the picture 
that arises in more formal derivations of the fluid / gravity 
correspondence in terms of derivative expansions of the 
field equations, where demanding normalizability at the 
boundary and regularity in the bulk (outside of black hole 
singularities) effectively constrains the gravitational dy- 
namics of the bulk to have as many "degrees of freedom" 
as the dual boundary fluid dynamics. For a review, see 
for example [15] . 

Applying the above reasoning to our metric fields g^ v , 
we construct regularized metric variables <? M „ that asymp- 
totically fall-off as ~ q: 



9tt 


= 9tt 4 


+ 


P)9tt 




9t P 


= 9t P ~\ 


-q 2 (l- 


+■ P) 2 9tp 




9t x 


= 9t x ~ 


1-9(1 + 


- P)9t x 




9pp 


= 9 PP - 


U(14 


- P)9pp 




9px 


= 9px~ 


fV(i 


+ p) 2 9px 




9xx 


= 9xx 


+-«(!- 


I" P)9xx 




996 


= 996 - 


1-9(1 -f 


P)(P 2 sin 2 x)9i> 




94>4> 


= 94>4> " 


f ff(H 


~P){p 2 sin 2 x sin 2 9)g n 


> (54) 



The term in the metric g^ that is conformal to S* 2 
can be kept track of by a single variable g^p , since we are 
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considering solutions that preserve an SO (3) symmetry 
that acts to rotate this S 2 



The extra factors p 2 sin 2 x 



and p sin \ sin 8 that appear in ( 54 ) are included to 



ensure regularity at the origin p — 0, for reasons similar 
to those discussed above. Using the regularized variables 
defined in ( 54 ) , the boundary conditions ( 14 ) can be fully 



captured by a simple set of Dirichlet boundary conditions 
at spatial infinity: 



9tt 


,=* 


= 


9t P 


P =i 


= 


9t x 


,=x 


= 


9 PP 




= 


9 PX 


p =l 


= 


9xx 




= 




,=* 


= 



(55) 



In polar-like coordinates, the Taylor expansion of ten- 
sorial quantities are typically either even or odd in p 
about p = 0; the extra factors of (1 + p) in ( [54] ) are to 
ensure that g^ and g pv have the same even/odd charac- 
ter in the limit p — > 0. We use standard results for the 
origin regularity conditions, which in our context read: 



Q p 9tt 


P - 


=0 


= 


d p9t x 


P - 


=0 


= 


dp9 P p 


P - 


=0 


= 


9 p9xx 


P - 


=0 


= 


d P g^ 


P - 


=0 


= 


9t P 


P - 


=0 


= 


9px 


P - 


=0 


= 



(56) 



Similar regularity conditions apply on axis: 



d x g tt 


X- 


= 0,7T 


= 


d x 9t P 


X 


= 0,TT 


= 


d x9 PP 


x= 


= 0,1T 


= 


®x9xx 


x= 


= 0,7T 


= 


9 X 9i: 


X 


= 0,7T 


= 


9t x 


x= 


= 0,TT 


= 


9px 


x= 


= 0,TT 


= 



(57) 



Elementary flatness imposes an additional relation 
among the metric variables on the axis, ensuring that no 
conical singularities arise there. Given our axial Killing 
vector d ( f >1 with norm-squared £ = g^, this statement is 
made precise by the condition [SD] 



4£ 



= 1, 



(58) 



which evaluates to the following in terms of our regular- 
ized metric variables: 



9XX\ X =0,7T ~ P W X =0,7T ' 



(59) 



We ensure that this relationship is satisfied by explicitly 
setting g xx in terms of g^p on the axis as dictated by ( 59 1 , 
instead of applying the regularity condition for g xx as 
displayed in (57). 



To find the appropriate regularized source function 
variables H^, we insert the above expressions (54 1 for 
the metric components into the definition ( 40 ) to find 



how the source functions deviate from their AdSs values 
near the boundary. Factoring out the appropriate 
powers of q so that H^{q = 0) = 0, and adding powers of 
(1 + p) to maintain the same expansions near the origin, 
we obtain 



H t = H t + q 3 {l + pfH t 
H p + q 2 (l + p) 2 H p 



H x + q 3 (l + pfH x 



(60) 



with boundary conditions: 



H t\ i 

L \ P =i 



H, 



p\ P =i 



H x\ p =i 

origin regularity conditions: 

d P^\ P=0 

d P H x\ p=0 

and axis regularity conditions: 



(61) 



(62) 



d x H t 


X 


=0,TT 


= 


d x H p 


X 


=0,7T 


= 


H x 


X 


=0,7T 


= 



(63) 



We also need a regularized massless scalar field variable 
(f> that asymptotically falls off as (j> ~ q, so we let 



= g 3 (l + p) 3 0. 
with boundary condition: 

origin regularity condition: 



(64) 
(65) 



d, 



X=0,7T 



(66) 
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and axis regularity condition: 



= ^(1) (f, X, 9, <t>)q + 4> {2) (t, X, 9, <j>)q 2 + 0(q 3 ) (70) 



d x4\ n = 0. 



C. Gauge Choice 



(67) 



Choosing generalized harmonic gauge conditions in 
AAdSs spacetimes is somewhat more subtle than in the 
usual asymptotically flat case. Roughly speaking, it 
turns out that it is not enough to simply demand that 
the metric and source functions satisfy the requisite rates 
of fall-off approaching the boundary as indicated in ( 54 ) 



and (60); rather, there are further restrictions amongst 



the leading-order behavior of certain fields that need to 
be explicitly enforced, so that the requisite fall-off is pre- 
served during evolution. 

To show this more clearly, we expand the regularized 
metric variables g pu , source functions H p , and scalar field 
<j) in power series about q = 

<W = 9(i)»»(t,X,0,<t>)<l+9(2)»v(t,X,9,(f>)q 2 +O(q 3 ) (68) 



and substitute these expressions into the field equa- 
tions (50). Since the GH form of the field equations re- 



sults in PDEs where the principle part of each equation is 
a wave operator acting on the metric (this fact guides our 
numerical solution method), we will schematically write 
this perturbative expansion of the field equations for the 
leading component ffnw of g^ v as follows 



//// 



(71) 



where we use the symbol □ to denote a wave-like opera- 
tor active within the (t, x, 9, <j>) subspace, and containing 
terms of the form cq ■d 2 g{\) P vldt l — C\ -d 2 g^i)^ v /dx 2 + 
Here, Co, C\, ... are coefficient functions that are in general 
different for each component of the field equations, but 
are regular and finite on the boundary. Their particular 
form is unimportant here, as we are interested in high- 
lighting the leading-order terms sourcing the wave-like 
equation on the RHS of (71). 



H P = H {1)fl (t, X , 9, </>)q + H (2)ll (t, x, 9, 



-0(q 3 ) (69) 



In this notation, the field equations near q = read 



□5(1)4* 

A<?(i)t P 

□S(i)*x 
□5(iW 
D 9(i)p X 

^5(i)xx 



_-2 



-0(q 

(-60p(i )tp - 

+ 2 9{i)t x ,x + 2 9{i) PP ,t ~ 5(i)xx, 
Oiq- 1 ) 

245(i) pp + 8<7 (1)xx 



(-8<?(i)p P + 4%) p )g 

3cotx5(i)t x + 245 (1)t - <7 ( i )tM 

;t - 2 5(l)V,t ~ 2 ^(l)p,t)? 



-2 



-<D(q 



( 



\i)tt 



(-QOg(i)px ~ 8 cot X5(i)xx + 8 cot X5(i)V + 24ff (i)x 
+9(i)tt, x ~ 2 9{i)t x ,t + 2 9(i)pp, x + 9{i)xx,x ~ 2 5(i)V:X ~ 2 %)p,x 



= (85(i) PP - 4iT ( i )p )<r 2 + Oiq- 1 ) 



165(i) + 16H (1}p )q- 2 + 0(q- 1 ) 

)q- 2 + 0(q- 1 ) 



= sin 



iin 2 x( 8 5(iW - ^H (1)p )q- 2 + O^ 1 ), 



(72) 



r 



and again we emphasize that □ is used to schematically 
denote a regular wave operator, though its specific form 
is in general different for each equation. For reference, we 



also list the leading-order behavior of the GH constraints 
(|45): 



Ct = {^9(i)t P - 8 #(i)t + 9(i)tt,t ~ 
C P = ( 8 5(i)« + 8 5(i) PP - 8 5(i)xx _ 
C x = {20g {1)px - 8H {1}x - 5(i)*t,x 



2 5(i)t x , x +9(i) PP} t +9{X)xx,t 
I6g (1)l p - 8H {1)p )q 3 + 0{q A ) 

+- 2 5(i)tx,t + 5(i) pp,x 9(i)xx,x 



2g {1 ^, t )q 4 + 0(q 5 ) 
a + 2g {1 ^, x )q 4 + 0(q 5 ) 



(73) 
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In order for the evolution to be consistent with each 



metric component's desired fall-off indicated in (54 1, the 
leading-order "s ource" on the right-hand side of each 



equation in (72) must scale as q° due to the Dirichlct 
boundary conditions (55 1. This implies that all terms 
of order q~ 2 and q^ 1 must vanish, and it would thus 
appear as though there were an additional hierarchy of 
constraints that need to be satisfied before the leading- 
order dynamics in the q° term becomes manifest (note 
that these constraints are not simply the GH constraints 



(73)). This is, in part, an artifact of having decomposed 
the field equations in a near-boundary expansion such as 
( [68] ): when solving the full equations consistently (i.e., 
with a Cauchy evolution scheme and initial data satisfy- 
ing the constraints, along with stable, consistent bound- 
ary conditions as discussed in the previous sections) , one 
usually expects that the evolution will "conspire" to pre- 
serve what appears as constraints in the perturbative ex- 
pansion. However, there are two potential complications 
in the AAdS case, highlighted in the above by the fact 
that the leading-order power in the expansions diverge as 
q-\ 

First, in the GH method, one is free to choose H p as 
the gauge, and though the structure of the field equations 
guarantee that the resultant solution will be consistent, 
there is no guarantee that a given choice of H p will pre- 
serve the desired asympotic fall-off of the metric (54). 



Case in point, suppose that we started with some initial 
value Hu(t — 0), and wanted to evolve to a gauge that 
in time becomes harmonic with respect to pure AdSs, 
namely H p (t > 0) = so that H M (t > 0) = H p 



in (60). Then (72) for g^tt immediately tells us that 
either 9m pp musTevolve to 0, or gmu ceases to remain 
regular. Whichever scenario unfolds, this choice of gauge 
leads to a representation of the metric that is not consis- 
tent with the fall-off assumed in (54), and as discussed 



in Sec. IV B it would generically be difficult (or even im- 
possible) to come up with a numerical scheme to stably 
evolve such a situation. 



Second, the form of the equations in (72) imply that 
regularity requires a delicate cancellation between terms 
in the near-boundary limit. Thus any truncation error 
introduced by a numerical scheme must be sufficiently 
small to effectively scale by some high power of q in the 
limit. In a typical finite difference scheme, the closest 
point to the boundary will be q = h , where h is the 
mesh-spacing there. Naively then, from (72) one would 



expect to need a discretization scheme that has local con- 
vergence order n + 2 to obtain global n th order conver- 
gence of the solution. Fortunately, this naive argument 
apparently does not hold with our code and the situations 
we have explored to date: we observe global second-order 
convergence using second-order accurate finite difference 
discretization. We are not sure why this is so; again, 
it could simply be that the near-boundary expansion is 
giving a misleading picture of the nature of evolution of 
the full set of equations, or it could be due to the par- 
ticular asymptotic gauge choice we use, described in the 



remainder of this section. 

We have experimented with a small handful of gauge 
choices which were unstable, including evolution to har- 
monic w.r.t AdS {H^{t > 0) — > 0, though as argued 
above one expects problems with this) , and a fixed gauge 
{dHu/dt — 0, and we note that our initial data is fully 
consistent with all the "constraints" in (72) and (73)). 



When a numerical solution is unstable, it is often dif- 
ficult to isolate the source: the gauge could be incon- 
sistent in the sense discussed above, the discretization 
scheme could be unstable for the particular set of equa- 
tions, there is a bug, etc. So here we simply list the 
asymptotic gauge condition we have empirically found to 
be stable (see Sec. VB for an explicit expression of the 



particular source functions used in this study), which in 
the notation introduced above is: 





5 




= 29{i)t P 


H (1)p 


= 2ff(i) pp 




5 




- 2 5(1)px 



(74) 



This choice was in part motivated by the asymptotic form 
of the field equations, in that these conditions, in conjuc- 



tion with the GH constraints (73) explicitly eliminate 
many of the order q~ 2 terms that appear in (72). How- 
ever, this is almost certainly not a unique choice for sta- 
bility, and one can anticipate that modifications would be 
required if, for example, the code is used to explore situa- 
tions with less symmetry, or to study scenarios where the 
boundary theory is "deformed" in a manner that alters 
the leading-order AAdS asymptotics. 



D. Initializing the 5-metric at t = 



Our goal in this section is to describe a choice of 
§tv\t=oi dt9tv\ t =o such that the source functions at t = 0, 



as evaluated via (40), are compatible with our target 
gauge (74). (The spatial components of the initial data 
<7ij| t=0 and dt9ij\ t= n come from the solution to the con- 



straint equations described in Sec. III). In the notation 
of the power-series expansion about q — performed in 



the previous section, the leading-order coefficients of ( 40 ) 
evaluated using (54) read 



5_ 1_ 1_ 1_ 

2#W*p + g9(i)tt.t - ^9(i)t x ,x + g9(i) PP ,t 

1_ 1_ 

+ gf(i)xx,t + l9(i)i>,t 

9(i)tt + 9{i) PP ~ 9(i)xx ~ 2 5(i)V 

5_ 1_ 1_ 1_ 

2~5(i)px ~~ g5(i)tt,x + 4?(i)tx,t + g9(i)p P ,x 

1 1 



\9{\)xx,x + l9{i)i>,x- 



(75) 



Inspection of ( 74 ) and ( 75 ) , together with our choice of 
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time symmetry (Qt5y| t _Q = 0), reveals that the following 
is required to leading order in the approach to p = 1 : 



9tt 1 1=0 = (Sp P + Sxx + 2 94> ) I t=o 
d t9tt\ t =o = 



(76) 
(77) 
(78) 



In practice we impose ( 77 1 and ( 78 ) over the entire com- 
putational domain, and impose ( 76 1 in an asymptotic 
sense: we set g~tt\ t=0 = for p < po, where po < 1 is a 



user-specified parameter, and smoothly transition to ( 76 ) 
as p — > 1. In this way, we keep the form of g tt\t = n as sim- 
ple as possible in the interior. The relations ( 74 1 and ( 75 ) 
leave the remaining metric variables gtp\ t=0 , dt9tp\ t=0 
and <?*xl t= o unconstrained, and we again take the sim- 
plest choice and set them to zero. In ADM language, 
these conditions amount to a choice of initial lapse, shift 
and their first time derivatives. 



V. NUMERICAL SCHEME 

The primary challenge in dealing with AAdS space- 
times lies in the nature of the boundary, which must 
largely be dealt with analytically, as discussed in detail in 
the previous section. Once this is done, the numerical dis- 
cretization and solution method is rather conventional. 
The methods we use are very similar to those described 
in |27j . and so we simply list them in Sec. VA referring 
the reader to [27] for the details. In Sec. VB we de- 
scribe the full (interior) form of the source functions that 
we use (in the previous section we only described their 



asymptotic behavior) , and in Sec. V C we show some con- 



vergence results from a representative solution. 



A. Numerical Solution Methods 

All equations and boundary conditions are discretized 
using second-order accurate finite difference methods, the 
only non-trivial boundary conditions being the Neumann 
conditions for the origin (ft = 0), and the axis (x = 0,7r) 
regularity conditions. The discretized Hamiltonian con- 
straint ( |36[ ) is solved using a full approximation storage 
(FAS) multigrid algorithm with v-cycling, and Newton- 
Gauss-Seidel iteration for the smoother. The evolution 
equations for the metric ( |50| and scalar field ( 15 ) (with 
energy- momentum tensor (34)), are discretized after sub- 



stituting in the definitions for our regularized metric vari- 
ables (54), scalar field (|64|), and source functions (60) 



The corresponding regular fields (g^,4>) are solved for 
using an iterative, Newton-Gauss-Seidel relaxation pro- 
cedure. In this study we have not explored any dynamical 
gauge evolution equations for H^, and simply set them 
to prescribed functions of the g^ as outlined in the next 
section. 

We use the excision method to solve for black hole 
spacetimes, which involves excising a portion of the grid 



within the apparent horizon (AH) , thus removing the ge- 
ometric singularity from the computational domain. Due 
to the causal structure of the spacetime within the hori- 
zon, all physical characteristics of the equations flow out 
of the domain (i.e. into the excised region). Thus exci- 
sion implies that the usual field equations are still solved 
on the excision surface, except that centered difference 
stencils are replaced with sideways stencils where appro- 
priate, in order to reference information that is only from 
the outside of the excision surface (in other words, no 
boundary conditions are placed there) . We search for the 
AH using a flow method; the excision surface is defined to 
be a surface with the same coordinate shape as the AH, 
but some fraction 1 — S ex smaller, where we typically use 
8 ex € [0.05,0.2]. One issue with using polar-like coordi- 
nates is that it incurs a rather severe restriction on the 
Courant-Friedrichs-Lewy (CFL) factor A that defines the 
time step At = A min(Ap, A%), where Ap and Ax are the 
mesh spacings in p and x respectively. Roughly speaking, 
with a uniform discretization in p and x, the condition 
for stability is A < p m in, where p m in is the smallest non- 
zero coordinate within the discrete domain. Ostensibly 
this occurs next to the origin, where p m in = Ap, so in 
the limit of high resolution the time-step size can become 
prohibitively small. For the tests and results presented 
here, we sidestep this issue by only studying black hole 
spacetimes, where excision removes the origin from the 
computational domain. 

Kreiss-Oliger dissipation [51] . with reduction of order 
approaching boundaries as described in [52], is used to 
help damp unphysical high-frequency solution compo- 
nents that sometimes arise at grid boundaries, in partic- 
ular the excision surface; we typically use a dissipation 
parameter of e = 0.35. 

We use the PAMR/AMRD 6 libraries to provide support 
for running in parallel on Linux computing clusters. The 
libraries also support adaptive mesh refinement (AMR), 
though all results described here are unigrid. 



Source Functions 



As described in Sec. IV C we are not free to arbitrarily 
choose the leading-order behavior of the source functions 
approaching the AdS boundary if subsequent evolution 
of the field equations is to preserve the desired asymp- 



totic form for the metric (14). Specifically, we want to 



consider a class of gauges that satisfy the asymptotic con- 
ditions (|74|). 



A naive implementation of ( 74 ) , wherein one would 
set H^lt > 0) everywhere on the grid via H t — 5/2g tp , 
H p = 2g pp , H x — 5/2<? px , would result in a discontinuous 
gauge at t — 0, as our method of solving for the initial 
data in general gives a different form for H^(t = 0) on 
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the interior. This is rather common in GH evolution 
i.e. the initial data provides source functions that are 
different from that of the target gauge, and the usual 
way to deal with this is to construct a hybrid gauge that 
smoothly transitions in time from the initial gauge to 
the target gauge. The specific transition we use is as 
follows. Denote the source functions coming from the 
initial data by = S^\ t _ , and define the functions 
F v to be the asymptotic constraints trivially extended 
into the interior 7 : 



Ft(t iP , X ) 
F p (t iP ,x) 
F x {t,p,x) 

Then, we choose 



^9t P (t,p,x) 

2g PP (t,p,x) 



with 



H M (t,p, X ) = H^(t,p, x )[g(t,p)} 

+ F v (t,p,x)[i-g(t,p)}, 



g(t,p) = exp(-z(i,p) 4 ), 



and 
f(p) 



z(t,p) 



1 



6/(p) + £i[i-/(p)]' 



1 — y 3 (6y 2 — 15y - 




10) 



P > 92 
P2> P> Pi 

otherwise 



(79) 



(80) 



(81) 



(82) 



(83) 



where y(p) = (p2~ p)/(P2~ Pi)), and £i,&,pi and p 2 are 
user-specified constants. The time-transition function g 
is such that g(0,p) = 1, g(0,p) = 0, g(0, p) — 0, and 
\im g(t,p) — 0; it is designed to give (80) the correct 



initial and target values, and transition between the two 
in a continuous fashion. The function /(p), on the other 
hand, is such that f(pi) = 0, f'(pi) — 0, f"(pi) = and 
f(p 2 ) = 0, f'(p 2 ) = 0, /"(pa) = 0; it is designed to let 
the transition occur with characteristic time £i for radii 
p < pi, interpolating to a characteristic time of £2 for 
radii p > p 2 . It is important near the boundary to reach 
the target gauge quickly, as this is where the delicate 
cancellations made possible by the gauge (74) are cru- 



cially needed. Accordingly, £ 2 is generally set to a small 
number. On the other hand, in may not be desirable to 
have such rapid gauge dynamics in the interior, and £1 
thus allows us to independently control the characteris- 



7 Of course, a more general implementation could allow for an 
arbitrary gauge in the interior, though for the spacetimes we 
have evolved to-date this simple choice has worked well. 



tic gauge evolution time there. On a typical run, we set 
pi = 0.0, p 2 = 0.95, & = 0.1, 6 = 0.0025. 



C. Convergence Tests 

To check the stability and consistency of our numeri- 
cal solutions we employ a pair of standard convergence 
tests. Here we show convergence results from one typi- 
cal representative case, namely strong scalar field initial 
data with non-trivial x-dependence. Specifically, the ini- 
tial data parameters (37) used were Aq = 10.0, Rq 
S 



0.0, 

0.2, w x — 4.0, w y = 16.0. The resulting evolu- 
tion describes a highly deformed black hole which settles 
down to an AdS-Schwarzschild solution with outermost 
apparent horizon radius rv, = 5.0. The physics of this 
solution will be discussed in Sec. |VI| together with addi- 
tional tests showing conservation of the boundary stress 
energy tensor in Sec. |VI C| 

First, in order to determine whether the evolution is 
stable and consistent, assuming that the solution admits 
a Richardson expansion we compute the rate of conver- 
gence Q(t, x l ) at each point on the grid for a given field 



Q(t, x l ) 



1 



ln(2) 



In 



Uh(t,x l ) - hh(t,x i ) 

f2h(t,X i )- /fc(t,S*) 



(84) 



Here, /a denotes one of g^v, H^, <f> from a simulation with 
mesh spacing A. Given that we use second-order accu- 
rate finite difference stencils, with 2 : 1 refinement in 
A between successive resolutions, and similarly in the 
time-step At, since we keep the CFL factor at a constant 
A = 0.2, we expect Q to asymptote to Q = 2 in the limit 
A -> 0. 

Second, to determine whether we are converging to the 
correct solution, i.e. to a solution of the Einstein field 
equations, we compute an independent residual of the 
field equations. This is obtained by taking the numerical 
solution and substituting it into a discretized version of 
G^u + A 5 g^ iv — 87rT M „. Since the numerical solution was 
found solving the GH form of the field equations, we do 
not expect the independent residual to be exactly zero; 
rather, if the solution is correct the independent residual 
should be purely numerical truncation error, and hence 
converge to zero. Thus, we can compute a convergence 
factor for it by using only two resolution results via 



QEFE(t,X l ) 



1 



ln(2) ™ V fh(t,*') 



In 



(85) 



Here, /a denotes a component of G^ v + A^g^ — 8ttT^. 
Again, given our second-order accurate finite difference 
stencils and with 2 : 1 refinement in A between successive 
resolutions, we expect Q to approach Q = 2 as A — > 0. 

Figs. [3] and [4] show L 2 -norms of (84) and (85) respec- 
tively, obtained from evolving the particular initial data 
described above. We used 4 different resolutions to help 
see the trends in the respective Q's. At early times, the 
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FIG. 3: Convergence factors (84 1 for the g pp grid function, 
constructed from a simulation run at 4 different resolutions; 
the highest resolution run has mesh spacing h/2. Here the 
L 2 -norm of the convergence factors are taken over the entire 
grid. The trends indicate that this grid function is converg- 
ing to second-order. Other grid functions exhibit similar 
trends. 
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FIG. 4: Convergence factors for the independent residual 
( 85 1 , constructed from simulations run at 4 different resolu- 
tions; the highest resolution run has mesh spacing h/2. At 
each point on the grid an L°° norm is taken over all compo- 
nents of the independent residual, and what is shown here 
is then the L 2 -norm of this over the entire grid. The trends 
in this plot indicate second-order convergence. 



solution is not yet in the asymptotic scaling regime for 
the convergence factors, however that they are typically 
greater than 2 indicates that the solution error is nev- 
ertheless small. Also, the trends going to higher resolu- 
tion, in particular at later times, appear consistent with 
second-order convergence. 

The early-time deviations in Fig. [4] coincide with an 
initial transient associated with the gauge transition de- 
scribed by ( 80 ) , that emanates from the boundary and 
dissipates as it travels further into the interior. 



This 

transient is clearly seen as a temporary dip in the point- 
wise convergence factors for independent residuals near 
the boundary; taking the L 2 -norm as was done in Fig. [5] 
to some extent masks this dip, though is still visible with 
the highest resolution curve in that plot. Nevertheless, 
this transient converges away in the sense that the region 
of the domain affected shrinks as resolution is increased. 



We now describe several early results obtained with 
this code. First, in Sec. |VI A[ we show solutions of the 
Hamiltonian constraint with a scalar field source ( 36 1 , 



demonstrating that this approach is capable of produc- 
ing initial data containing trapped surfaces. Then, in 
Sec. IVIBI we show the evolution of initial data describ- 
ing highly deformed black holes that subsequently shed 
their asymmetries via quasi-normal ringdown. A proper 
extraction of quasi-normal modes, and in particular a 
meaningful comparison with perturbative calculations 
where these modes can be defined, requires the identi- 
fication of a reference background metric and a transfor- 
mation of the solution to a gauge consistent with that 
of the perturbative calculations. We have made no at- 
tempt in this direction, besides matching the areal ra- 
dius on the extraction sphere, and do in fact see what 
appear to be gauge modes. Nevertheless, we can extract 
the leading order linear quasi-normal modes, and for the 
higher angular number modes, use a simple forced har- 
monic oscillator model to identify what appear to be non- 
linearly excited harmonics of the lower angular modes. In 
Sec. |VI C] we discuss the extracted boundary stress energy 
tensor of these solutions, and in Sec. |VID| analyse how 
well the CFT state can be described by hydrodynamics. 
We find that the extracted stress energy tensor is con- 
sistent, essentially to within numerical truncation error, 
with that of a viscous, conformal fluid from t = on- 
wards. In Sec. |VIF| we transform these solutions onto 
a Minkowski piece of the boundary, and find an initial 
fluid geometry that resembles a Lorentz-flattcncd pan- 
cake. Defining the beam-line direction as the one along 
which the initial data is flattened, the evolution of this 
fluid exhibits both longitudinal and transverse flow rel- 
ative to the beam-line, with most of the energy flowing 
along the longitudinal directions. 



Strong-field Solutions of the Hamiltonian 
Constraint 



To generate black holes spacetimes, we choose initial 
data where the deviation from pure AdS5 is sourced by 
a highly compact distribution of scalar field energy, with 
profile given by (37). In this paper we only consider 
free, massless scalar fields, so V((j>) — in (35). In this 
section, we begin by focusing on spherically symmetric 
initial data, so w x = w y = 1 in (37). 

Fig. [5| su mmarizes solutions of the Hamiltonian con- 
straint ( 36 ) by plotting the maximum value of the confor- 



mal factor versus the maximum of the scalar field (both 
maxima occur at the origin of the domain for these initial 
data). Let us begin by contrasting the qualitative behav- 
ior of this plot with its counterpart in the asymptotically 
flat case, presented in [53] - This earlier work employed 
the conformal thin sandwich method to solve the con- 
straints. There, it was discovered that there exists a crit- 
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ical point above which no numerical solutions for the ini- 
tial data could be found, and where the conformal factor 
diverges as the amplitude approaches the critical point. 
Even more surprisingly, it was found that generalizing to 
the extended conformal thin sandwich method gives rise 
to a branch point instead of a critical point, with solu- 
tions along the upper branch exhibiting non-uniqueness 
for any given set of initial data. 
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FIG. 5: Maximum conformal factor vs. maximum ampli- 
tude of initial scalar matter distribution, with Gaussian pro- 
file (|35} with A = max , 8 = 0.2 (with AdS scale L = 1), 
w x — w y = 1 and Ro = 0. We differentiate between 
"strong-field" and "weak-field" data based on whether there 
is a trapped surface present on the initial slice or not, re- 
spectively. (Though of course this distinction is somewhat 
arbitrary, particularly since subsequent evolution of weak- 
field data could eventually result in black hole formation, 
as argued in [19] even for arbitrarily small amplitude initial 
data.) The value of 4>max beyond which trapped surfaces are 
found in the initial data is indicated by the dashed vertical 
line. The open circles denote numerical solutions, while the 
solid lines are fits to the data, as shown. 
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FIG. 6: AdS mass vs. maximum amplitude of massless 
scalar matter, with Gaussian profile defined by A — (i maj! , 
(xq, j/o) = (0, 0), 8 = 0.2. Conventions are the same as those 
used in Fig. [5] In cases where trapped surfaces are present 
(to the right of the dashed vertical line), an estimate of 
the mass based on the area of the apparent horizon gives a 
value very close the asymptotic mass plotted here, though 
systematically smaller (the two values would essentially be 
indistinguishable on the scale of this figure, hence we only 
show the asymptotic mass to avoid clutter). 



uniqueness of solutions to non- linear elliptic PDEs can be 
understood by applying a maximum principle, where the 
signs and relative magnitudes of coefficients in the PDE 
are crucial; one can show that in this regard cosmological 
the constant in AAdS spacetimes "helps" . 

Fig. [6] shows a plot of the conserved mass M of the 
spacetime versus scalar field amplitude, computed from 
the quasi-local stress tensor (21) as follows. We take 
a spatial t = const, (here t = 0) slice in dM q , with 



induced 3-metric o"u„, lapse N and shift N l such that 
T.^dx^dx" = -N 2 dt 2 + a ij {dx i + N l dt)(dx j + NHt), 
then compute 



In contrast, we find no divergent behavior in the 
conformal factor, nor any behavior suggesting the non- 
uniqueness of solutions 4> m ax , at least in the regime where 
the cosmological length scale is relevant 8 . Furthermore, 
the linear dependence of the maximum of the confor- 
mal factor versus amplitude for large amplitude data 
suggests that the AdSs Hamiltonian constraint admits 
conformal solutions for arbitrarily strong initial matter 
distributions. Of course, given the presence of the cos- 
mological constant and its relevance on these scales, it 
is not too surprising that we find such qualitatively dif- 
ferent behavior compared to the asymptotically flat case. 
From a more formal perspective, the local existence and 



We have not investigated the limit where the characteristic scale 
in the initial data <5 is much less than the cosmological scale L, 
where one might intuitively expect L to become irrelevant in 
governing the local nature of the solution, and results consistent 
with the asymptotically flat case 1531 may be recovered. 



M = lim 



(86) 



where u M is the time-like unit vector normal to t = const. 
Note that for a vacuum AdS-Schwarzschild black hole, 
this prescription gives a result consistent with the usual 
definition of its mass from the analytic solution, namely 
M = 37r(r 2 /8), where r = r#\/l + r 2 H /L 2 , given a 
horizon with areal radius ru- 



B. Quasinormal Ringing 

As a first application of our evolution scheme, we study 
the quasi-normal ringdown of an initially high distorted 
(i.e. non-spherical) black hole. Here we focus on the 
dynamics of the bulk, and in the following sections dis- 
cuss the corresponding dynamics of the CFT boundary 
stress tensor. For some previous work on the subject 
see for example [351 [33J I53H57] . For a review of black 
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hole quasi-normal modes in AdS see [55]. Again, we use 
large-amplitude scalar field data to create an initial slice 
of the spacetime containing a trapped surface, and here 
introduce a non-trivial x-dependence by adjusting the 
shape of the profile through the parameters w x and w y 
as defined in (37 1. Specifically, unless otherwise stated 
we choose i? = 0.0, S — 0.2, w x — 4.0, w y — 16.0, and 
vary the amplitude Aq to control the size of the resultant 
black hole. 

First, to demonstrate that we are looking at relatively 
large perturbations of the spherical black hole, in Fig. [7] 
we plot the ratio of the equatorial c eq to polar c p proper 
circumferences of the apparent horizon versus time 9 . A 
geometric sphere has c eq /c p = 1, whereas a geometric 
disk c eq /cp = it/2. The w y /w x — 4 case studied in de- 
tail here has c eq /c p (t = 0) w 1.2, so a fairly sizeable 
deformation. For certain quasi-normal mode and hy- 
drodynamic extraction results we also show data from 
the Wy/w x = 32 case, which has a very large initial 
c eq /c p (t = 0) « 2.2. 
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FIG. 7: The ratio c eq /c p of proper equatorial to polar cir- 
cumference of the apparent horizon versus time. w y /w x 
denotes parameters in the initial data describing the scalar 
field profile (37 1, and w y /w x = 4 is the canonical case stud- 



ied in this paper, though for interest we also show examples 
representing a weaker and stronger initial asymmetry. In 
all cases, the final black hole has radius th ~ 5. Note that a 
geometric sphere has c eq /c p — 1, whereas a geometric disk 
has Ceg/cp = it/2. Thus, even the w y /w x = 4 case is ini- 
tially a rather large deformation of the 3-sphere; the curve 
for the more extreme case w y /w x = 32 implies that at early 
times the Gaussian curvature is negative over at least some 
regions of the horizon. 



9 Note that the apparent horizon is to some extent slicing depen- 
dent. However, given the symmetries in our problem, in partic- 
ular that t = is a moment of time-symmetry (where we see the 
largest deformation of the horizon), it is difficult to imagine how 
the intrinsic geometry of the apparent horizon does not give a 
good indication of the relative magnitude of the spacetime per- 
turbation. Nevertheless, it would be interesting to compare to 
the event horizon, though finding event horizons are more com- 
plicated and we leave that to a future study. 



To extract the quasi-normal modes, the first step 
would be to transform our metric to a gauge consistent 
with standard perturbative calculations for quasi-normal 
modes. This is a rather non-trivial step in general, and 
we do not do so here, hoping that our asymptotic gauge 
choice is sufficiently close to that of the perturbative cal- 
culations that artifacts introduced are small. As we shall 
see, this seems to hold to good approximation, though 
we do find modes that appear to be pure gauge. 

We begin by projecting our field variables onto the 
spherical harmonics on S 3 ; see [59] 
standard spherical harmonics on S 2 



We first need the 



121 + 1 {l-m)\ 
4tt (7 + to)! 



P ; m (cos6»)exp(im(/)) (87) 



where the associated Legendre polynomials are defined 
as: 

m ( — I)™ S l + m ^ 2 ; 

The scalar spherical harmonics on S 3 are then: 



S(klm) = (-1)¥(2Z)!! 



'2(* + l)(fc-0 



7r(fc-M + l)! 
xC l + 1 l (cosx)sm l (x)Y lm (e, ( f>) (89) 

where the Gegenbauer polynomials are defined as: 



C a k {x) = 



(-2) fc r(fc + q)r(fc + 2 Q ) 2] - a+ l/2 
k\ T(a)T(2k + 2a) [ ' 



xf^(l-* 2 ) fe+a - 1/2 - 
ax K 



(90) 



Our restriction to solutions that preserve a symmetry in 
9, 4> mandates that we have / = 0, m = 0. 

For perturbations of the scalar field, these S(klm) 
scalar harmonics exist for all k > 0. Fig. [8] shows the 
scalar field projected onto the 3-sphere (k = 0), for a 
simulation with initial w y /w x = 4 and whose final state 
black hole has horizon radius = 5. The quasi-normal 
mode frequencies extracted from the scalar field are col- 
lected in Appendix [C] under Tabl e |IV| for the n = 
fundamental mode and under Table IVTfor the n — 1 first 
overtone. This was done for several r% cases, wherein 
each pair of fundamental and first overtone was found by 
a simultaneous non-linear least-squares fit to two damped 
sinusoids, one corresponding to the n = fundamental 
and the other corresponding to the n — 1 overtone. Each 
sinusoid is described by four parameters 10 . The n = 
fundamental frequencies uj r and scale linearly with r^, 



10 These parameters are amplitude A, phase ip, imaginary frequency 
u)i, and real frequency u> r , so the fitting functions take the form 
A exp(— LJit) cos(u r t + tp). 
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and their dependence on k shows that u> r increases with 
k and that u>i decreases with k. The n = 1 first overtone 
frequencies uj r and Wj also scale linearly with r/,, as with 
the n = frequencies; however, the dependence on k dif- 
fers from the n = fundamental case in that the uj r and 
Wi both increase with k. To make contact with earlier 
work on quasi-normal modes in AdS, we note that the 
frequencies extracted for the n = fundamental are a 
close match to those found in the seminal study |54| on 
this subject. 
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FIG. 8: The leading-order behavior of the scalar field vari- 
able <f> near the boundary 5 = 1, projected onto the 3-sphere 
by $ fe=0 = / dfi 3 (0/g)S(OOO), and plotted over a global time 
interval of t S [0, 0.62] (open circles). A fit (solid line) 3>{!fo 
is extracted using the data inbetween the dashed vertical 
lines. The inset shows a logarithmic plot of the data and 
the fit over the full global time interval. 
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FIG. 9: The leading-order behavior of the metric variable 
g xx near the boundary q = 0, projected onto the §(200) 
scalar harmonic by (Q xx ) k = 2 = / d^l^(<j xx l 'g)§(200), plot- 
ted over a global time interval of t £ [0, 87r] (open circles). 
A fit (solid line) {Q 

xx)/d2 I s extracted using the data in- 
between the dashed vertical lines. The inset shows a loga- 
rithmic plot of the data and the fit over the full global time 
interval. Other metric variables show similar behavior. 
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TABLE I: The fundamental (n = 0) quasi-normal mode 
frequencies u r — iui extracted from the metric variable g xx . 
These are shown for various horizon radii rh of the final 
state AdS-Schwarzschild black hole and for 5*0(4) quantum 
numbers k — 2, 1 = 0, m — 0. Uncertainties are estimated 
from convergence studies. 
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TABLE II: The first overtone (n = 1) quasi-normal mode 
frequencies u r — iuii extracted from the metric variable g xx . 
These are shown for various horizon radii rh of the final 
state AdS-Schwarzschild black hole and for 5*0(4) quantum 
numbers k — 2, 1 = 0, m — 0. Uncertainties are estimated 
from convergence studies. 

For the metric, the S(klm) scalar harmonics are re- 
stricted by < I < k. The k = 1 mode is not associated 
with any physical degrees of freedom |60j . and these can 
be ignored; indeed our choice of initial data where the 
scalar field profile is symmetric about \ = 7r/2 prevent 
these from being excited 11 . The k = mode corresponds 
to a perturbation of a black hole that is itself spherically- 
symmetric in all the S 3 angles %, 6, <j>. Thus, as a con- 
sequence of Birkhoff 's theorem, we might also expect to 



In addition to the scalar harmonics, the metric also in general 
admits vector Y(klm) (1 < I < fc) and tensor harmonics T(fcZm) 
(2 < £ < k), but our 50(3) symmetry precludes the excitation of 
these modes: the vector and tensor harmonics require 8 depen- 
dence (i.e. I ^ 0) that our symmetry does not allow. 
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ignore this k — case as we have with k = 1. How- 
ever, despite this reasoning, we see non-trivial dynamics 
in the data's projection onto the k = scalar harmonic. 
To understand why, recall that we have not transformed 
our metric to a gauge that is consistent with the stan- 
dard one assumed by the pcrturbativc calculations for 
quasi-normal modes. Consequently, the metric compo- 
nents that we extract generally contain both physical 
quasi-normal modes as well as contributions from gauge 
modes that are introduced by our "non-standard" choice 
of time-slicing, and of the r and x coordinates on each 
slice. This gauge contribution is of course present in the 
data for all fc, and is most obvious in the k — pro- 
jection since it contains no physical quasi-normal modes. 
For k > 2, this gauge contribution manifests itself as a 
decaying mode (i.e. a mode whose frequency has a neg- 
ligibly small real part). In the following, we consider the 
k > 2 modes of the metric, taking into account this gauge 
contribution. 

Fig. [9] shows a representative metric variable projected 
onto the §(200) (fc = 2) scalar harmonic, for a simulation 
with final state horizon radius r% = h and initial asym- 
metry w y /w x — 4. The quasi-normal mode frequencies 
extracted from this representative metric variable are dis- 
played in Table H] This analysis was repeated for several 
simulations with varying but fixed w v /w x = 4. Fits 
to damped sinusoids yield n = fundamental frequencies 
with imaginary parts u>i that scale as ~ l/^/i, and real 
parts ui r that are largely insensitive to these changes in 
rh. 

Obtaining the §(200) harmonic's n = 1 overtone from 
the metric variables is more difficult, largely because it 
decays much more quickly than the n = fundamental. 
To obtain good fits, we focus on an early-time segment 
of the data, and substract the n — fit, then fit to the 
remainder. The quasi-normal frequencies extracted in 
this way are tabulated in Table |Tlj This table shows that 
the n = 1 fundamental frequencies ui r and uii both scale 
linearly with rh- To compare with earlier work, notice 
that these extracted frequencies for the n = 1 overtone 
are a close match to the first set of "fast-modes" found 
in |33j . and that the extracted frequencies for the n = 
fundamental closely match the low-lying "slow-modes" 
found in the same study. 

Discrepancies with the linear quasi-normal mode de- 
scription start to appear in the next highest scalar har- 
monic §(400) (A; = 4): direct fitting yields a dominant 
frequency that does not match 33! ■ The k = 4 fun- 
damental mode with frequency w 4 ps 2.949 + i3.428/r^ 
is present, but is overshadowed by a mode with fre- 
quency cjf 1 fa 3.312 + £1.627/77,, which is close to double 
that of the k = 2 fundamental UJ2 ~ 1.652 + «0.826/r^ 
that appears in Table [TJ We expect that this frequency- 
doubling in the k = 4 modes arises from a non-linear 
mode-coupling, which we attempt to model as a damped 
harmonic oscillator driven at double the frequency of the 



k = 2 fundamental mode 12 . Given a k = 2 fundamental 
mode ^2(0 = A2 exp (— iw^t), the k = 4 mode \&4 (t) in 
this simple model satisfies 

<9 t 2 * 4 + kj^ 4 - A 4 a t * 4 = B exp (-2iu 2 t) , (91) 



A4 = tJa, k\ = (u>4 r ) + (u>a) (with the r and i sub- 
scripts denoting the real and imaginary components of 
the corresponding number respectively), and we expect 
the driving amplitude B to scale as B ~ (A2) 2 . 

4 fundamental 



The solution of ( 91 ) is a sum of the k 



mode and the driven mode 



* 4 (t) = A± exp (-iuit) + Af l exp (-2£cj 2 *) , (92) 

where A± is a constant depending on the initial data, and 
Af bl is a (complex) constant depending upon the other 
parameters of the model via 



Af l = {Q r +iQi)- L B, 



(93) 



where we have introduced for notational convenience 
Wi = 4a; 2r (2cLi 24 - u)u) and Q r = (w 4r ) 2 + (w 4l ) 2 - 

4 ((W 2 r) 2 - (^2i) 2 + W4jW2i). 

Our strategy is to compare the predictions of this sim- 
plified model of non-linear mode-mixing to the parame- 
ters extracted from the fits. In particular ((93|) gives us 



quantitative predictions for the relative amplitude and 
phase difference between the frequency-doubled k = 4 
mode and the k = 2 fundamental mode that sources it: 
that the phase difference should be A<p = arctan(w i /w r .), 
and that Af bl should scale as (A2) 2 . To test these predic- 
tions, we fit the k = 4 data to two damped sinusoids si- 
multaneously (in addition to the decaying gauge mode), 
though fixing their frequencies to u; 4 and 2w2- We can 
then check if the extracted amplitude and phase of the 
frequency doubled mode matches the model prediction. 
The results, presented in Tabic |III| for three values of 
Wy/w x , and an example fit is shown in Fig. [10| for the 
Wy/w x = 4 case, show good agreement with the model. 

To quantify precisely how well the fits describe the 
actual data, we compute a residual that measures the 
difference between the data and the fits. This residual 
quantifies the part of the dynamics that we have not 
been able to fit to with linear quasi-normal modes sup- 
plemented by the pure decay gauge mode, and in the 
k = 4 case, the mode arising from non-linear mode- 
mixing modeled by (91 1. The solid blue curve of the 
first panel of Fig. [TT] depicts a normalized residual that 
is generated from the fit shown in Fig. [9] To illustrate 
the dependence on w y /w x , this analysis is repeated for 
the largest and smallest w y /w x we have considered. The 
solid blue curve shows that the fit to the k = 2 fundamen- 
tal quasi-normal mode for t > 1 captures all but ~ 1% 



12 This is reasonable if one considers the form the field equations 
take, perturbing about a black hole solution to second-order. 
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Adbl 
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U] 9 






l.i 

4 

32 


(1.49 ±0.09) x 10~ 3 

2.75 ±0.01 

(5.39 ±0.04) x 10 1 


(6.44 ±0.02) x 10~ 3 

9.62 ±0.01 

(1.16 ±0.01) x 10 2 


(4.40 ±0.03) x 10~ 4 
0.508 ±0.030 
1.59 ±0.06 


(1.23 ±0.23)^ 
(1.26 ±0.04)^ 
(1.16 ±0.01)^ 


0.832 ±0.011 
0.873 ±0.032 
0.894 ±0.066 


4.03 ±0.02 
4.00 ±0.02 
4.60 ±0.02 



TABLE III: Extracted parameters from the k — 4 fit corresponding to Fig. [lOjat w y /w x — 4, as well as for the 
largest and smallest w y /w x cases considered. The fit includes the fundamental k = 4 mode with amplitude A±, 
a mode of amplitude Af bl that has twice the frequency of the k = 2 fundamental mode of amplitude A2 , and a 
putative gauge mode A\ exp(— wf^t). The model for the mode coupling discussed in the text predicts a phase 
difference between the frequency-doubled k — 4 mode and its k — 2 fundamental mode source of Aip ~ 0.842, 
and that A4 scales as the square A2 ; these two quantities are shown in the last two columns of the table, and 
indicate the solution is reasonably consistent with the model. Uncertainties are estimated via convergence. 
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FIG. 10: The leading-order behavior of the metric variable 
g xx near the boundary q = 0, projected onto the §(400) 
scalar harmonic by (Q xx ) k=4 = J dQ3(g xx /q)§(400), plot- 
ted over a global time interval of t £ [0, 87r] (open circles). 
A fit (solid line) (Q 

xx)fcL4 i s extracted using the data inbe- 
tween the dashed vertical lines. This is a fit to a gauge 
mode, and two damped sinusoids with fixed frequencies 
loa « 2.949 + i3A28/r h and 2 uj 2 ~ 2 x (1.652 + i0.826/r h ) 
(see Table |III| for the corresponding parameters of the fit). 
The inset shows a logarithmic plot of the data and the fit 
over the full global time interval. Other metric variables 
show similar behavior. 



of the perturbation projected onto the §(200) harmonic, 
for the w y /w x — 4 case. The other curves show that a 
similar fit does worse for larger w y /w x , which is consis- 
tent with the expectation that other unmodeled effects 
(i.e. those not accounted for by the linear quasi-normal 
modes and our simplified model of the gauge and non- 
linear effects) should become more significant for data 
with larger deformations. Similarly, the solid blue curve 
of the second panel in Fig. [TT] depicts the normalized 
residual corresponding to the fit shown in Fig. [TUJ It 
shows that for t > 1, the fit to the fixed- frequency k = 4 
fundamental mode, the frequency-doubled k — 2 mode 
and single gauge mode captures all but ps 5% of the solu- 
tion projected onto the §(400) harmonic for w y /w x — 4. 
Again, the corresponding residual grows with larger ini- 



tial asymmetry. 

Finally, to quantify how much of the full solution is not 
accounted for by the §(200) and §(400) harmonics that 
we discussed above, we look at the normalized difference 
between the square of the full solution integrated over the 
3-sphere, and the sum of the squares of all projections 
onto the §(fc00) harmonics for all k < k max - The result 
of this procedure is summarized in Fig. [12] 



C. Boundary Stress Tensor 



In terms of the regularized metric variables ( 54 ) , and 



setting L = 1, the stress energy tensor of the boundary 
CFT (Il9j) evaluates to 



( T tt) cft 

CFT 
(Txx) CFT 

(Tee) cft 



64tt 
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2^ 
1 



64tt 

. 2 
sin \ 



(-24g pPtP - 32g xx ^ p - 64^) 



i-9t x , P ) 



(-32g tt „ + 2Ag i 



pp,p 



64 6V>„ 



64tt 



(-32g U:P + 24g pp , p + 32& 



(94) 



CFT 



sin 9 (Tgg) CFT . In the remainder of this 



and (T00) 

section we discuss the CFT stress tensor extracted from 
a representative numerical solution, namely the r% = 5 
quasi-normal ringdown spacetime described in the previ- 
ous section. 



First, for a consistency check of expressions (94), we 



test whether (to within truncation error) the stress tensor 
is traceless and whether it is conserved with respect to 
the Levi-Civita connection on the Rx5 3 boundary. Re- 
sults for the trace and two non-trivial components of the 



divergence are displayed in Fig. 13 and 14 respectively. 
These plots demonstrate that as resolution is increased, 
we are indeed converging to a CFT stress tensor that is 
conserved and traceless, i.e. to matter that obeys the hy- 
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FIG. 11: The normalized residuals for the k — 2 and k — 4 
fits, corresponding to Fig. [9] and Fig. 1 10| respectively, where 

(K xx ) k = {{Qxx)k - (Sxx)fc 1 *) / |(5xx)fc l ')|- This quantity 
shows the extent to which the projection of the metric el- 
ement g xx onto the S(fc00) harmonic fails to be described 
by the sum of the k quasi-normal mode, gauge mode, and 
for the k — 4 case the non-linear mode. At times before 
( « 1, the residual becomes large going to smaller t due 
to a growing phase offset between the data and the fit (see 
Fig. [9] for the k = 2 case, and Fig. 10 for the k — 4 case), 



and at very early times there is an additional contribution 
from the early-time transient. At late times there is a slight 
increase in the residual. This is largely a consequence of 
the fitting procedure and that we are plotting a normalized 
residual : the fit is dominated by the data at early times 
as the perturbation is decaying exponentially, hence a small 
phase difference between the fit and data will result in a 
normalized residual that grows with time. 



FIG. 12: The normalized difference (T- , xx)k ma x = 
(/ dn 3 (g xx /qf -YXSr^Qxxfu) 1 1/ dn 3 (g xx /q) 2 \. For 
the Wy/tVx = 4 case, the contribution due to the k > 2 
(k > 4) modes constitutes w %2 (%0.1) of the metric per- 
turbation, and decreases with time (the metric modes decay 
faster with higher k). 




drodynamic equations of motion, and whose equation of 
state is consistent with conformal invariance. Note that 
the early-time transient related to the initial gauge dy- 



namics discussed in Sec. V C is also visible in these plots, 
though as with the independent residual it also converges 
away in the sense that it occupies a smaller region of the 
spacetime domain with increased resolution. 



In Figs. 15 and 16 we display representative compo- 



nents of the boundary CFT stress tensor on K x S space- 
time diagrams, specifically the energy density (T tt ) CFT 
and S* 2 component of the pressure (Tgg) CFT respectively. 



FIG. 13: The trace (r'V)cFT °f the boundary stress ten- 
sor, constructed from a simulation run at 4 different reso- 
lutions, labeled by the mesh spacing relative to the highest 
resolution run h. Here the L 2 -norm of (T' A ' M } CPT is taken 
over the entire grid, and the trends indicate convergence to 
a trace-free stress tensor with increasing resolution. 



At t = the boundary state is clearly inhomogeneous, 
and in time evolves to a homogeneous state in a manner 
that mirrors the quasi-normal decay of the spacetime to 
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FIG. 14: Two components of the divergence of the stress 
tensor, V M <T Mt ) 

cft (t°p) an d V (Tmx)gft (bottom), con- 
structed from a simulation run at 4 different resolutions, 
labeled by the mesh spacing relative to the highest resolu- 
tion run h/2. The L 2 -norm of the respective components 
are taken over the entire grid, and the trends indicate con- 
vergence to a conserved stress tensor with increasing reso- 
lution. 



a static black hole. 



D. Hydrodynamic Description of the Boundary 
CFT 



Essentially all known physical systems at sufficiently 
high temperature, including those described by quantum 
field theories, exhibit hydrodynamic behavior once local 
thermodynamic equilibrium has been attained. In the 
gauge/gravity duality, stationary black holes are dual to 
equilibrium thermal states, and studies have shown that 
perturbations of the bulk spacetime manifest as hydro- 
dynamic behavior in the boundary CFT. See |49j for a 
recent review, and [61 for a review of the "blackfolds" 
approach, which is similar in spirit but connects the dy- 
namics of a perturbed black brane with that of an en- 
closing world-volume via a derivative expansion of the 
Einstein field equations. Here, we have studied the evo- 
lution of initially highly distorted black holes, far from 
the perturbative regime, and so the question naturally 
arises: at what time does hydrodynamics become a good 
description of the boundary state? More precisely, we 
would like to know whether the extracted (T) ll/ ) CFT on 
the Rx5 3 boundary is that of the particular kind of fluid 
predicted by the duality, namely an J\f — 4 SYM fluid 
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FIG. 15: The energy density {T tt ) CFT of the boundary CFT, 
displayed on a (t, x) spacetime diagram, extracted from the 
Th = 5 quasi-normal black hole ringdown simulation de- 
scribed in Sec. |VI B| The first panel corresponds to the 
Wy/w x = 4 case, and the w y /w x — 32 case in the second 
panel is included for comparison. 



with equation of state e = 3-P (where e and P are the 
energy density and isotropic pressure in the rest frame 
of the fluid, respectively), and transport coefficients that 
match those found via holographic methods [6"2"] . 

In the previous section, we showed that convergence 
trends in the solution imply that in the continuum limit, 
(Tjw) CFT is conserved and traceless. Thus, two neces- 
sary ingredients are already satisfied: conservation is re- 
quired for the effective hydrodynamic variables to satisfy 
the Navier-Stokes equations, and for an isotropic fluid, 
tracelessness implies e = 3P. In this section we show 
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FIG. 16: The S 2 component {T ee ) CFT of the boundary CFT 
stress tensor, displayed on a (t, \) spacetime diagram, ex- 
tracted from the = 5 quasi-normal black hole ringdown 
simulation described in Sec. |VI B| The first panel corre- 
sponds to the Wy/w x =4 case, and the w y /w x = 32 case in 
the second panel is included for comparison. 



that one can map (T^ i/ ) CFT to a corresponding set of 
hydrodynamic variables, and that essentially from the 
initial time, their behavior is consistent with that of an 
Af = 4 SYM fluid, at least to within truncation error. 

In a low-energy effective description, the hydrody- 
namic stress tensor can be expressed as a velocity- 
gradient expansion, i.e. as a power series in the covariant 
gradients of a local fluid 4- velocity one form u v : 



Tftv = (e + -P)u M u„ + Pg^v - 



n 



(95) 



shear tensor 



a „u =±lia±v p V{aUfj) _ — V a u a U*", (96) 

and subsumed all higher-order terms under II Mt/ . All co- 
variant differentiation is performed with respect to the 
boundary metric g^dx^dx" = —dt 2 + d\ 2 + sin 2 x(d9 2 + 
sin 2 8dcf> 2 ), and _L MI/ in (96) is the projector onto the spa- 



tial slices orthogonal to the fluid 4-velocity 



_L' 



ILV | Li V 

+ wu . 



(97) 



In this sub-section, we will for now ignore the higher- 
order terms, i.e. we set = 0. In terms of mapping 
from stress tensor to hydrodynamic variables this should 
be a good approximation except in situations where the 
shear cr M „ of the flow becomes small, which does occur 
periodically during the evolution of the distorted black 
holes discussed here. These higher order terms would 
lead to additional problems regarding the uniqueness of 
the mapping that we will now attempt. In the sub-section 
following this one, we will employ a different strategy and 
will instead extract a minimal subset, then, assuming 
an Af = 4 SYM fluid, test for consistency with higher 
order transport coefficients (instead of trying to directly 
measure all these quantities, as we will now proceed to 
do). 

Let us first identify a set of independent fluid vari- 
ables which we can use to characterize the stress ten- 
sor. This set certainly includes e and P, the energy 
density and isotropic pressure in the rest frame of the 
fluid, respectively. Given the SO (3) symmetry of our 
solutions, the (unit) 4-velocity vector must take the 
form = 7(1, v, 0,0) in (t, %, 9, <fi) coordinates, with 
7 = 1/vi — v 2 . This gives us a third fluid variable to 
add to our list: v, the x-coordinate velocity of the flow. 
As for the shear a^, notice that it is symmetric, trace- 
less and satisfies the identity u^a^ u — 0. Together with 
the imposed 5*0(3) symmetry, these imply that a „ v only 
has one degree of freedom in d = 4 dimensions. Thus 
defining a xx = a, one can straightforwardly show that 
the only other non-zero components of cr^„ are: 



-vcr, <Jtt 



v 2 a, 



a ee — . 7 
sin 1 



sin 2 x 



2 7 2 



(98) 



where we have introduced the (symmetric, traceless) 



Given how 77 and a always appear as a product in (J95J), 
we treat the two as a single quantity rja for the purposes 
of extraction. Thus, ignoring the higher-order terms II M „, 
the variables (e, P, v, 77a), each of which is a function of 
(i, x) m general, completely describes a conformal fluid 
flow on M x 5 3 that preserves our SO (3) symmetry. 

On the gravity side, the quantities we measure are the 
components of the boundary stress tensor, which we de- 
note T^v = (7) t j / ) CFT henceforth for the sake of brevity. 
In SO(3) symmetry, there are 5 non-zero stress tensor 
components, 4 of which are independent. We can relate 
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these 4 to the fluid variables via ( 95 1 



T tt = (e + P) 
Tt x = ~{e + P) 



1 



T, 



xx 

Too 
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1-v 2 

V 

1-v 2 
v 2 



- P - 2-qav 1 
2r/av 



P - 2r/er 



1 -v 2 
sin 2 x (P + ?/cr(l - v 2 )) 



(99) 



where — sin 8Tgg. Inverting these relations in favor 
of the rest-frame hydrodynamic quantities (e, P, v, rye), 
and defining the auxiliary quantities 



T$ = 
we obtain 



£ = \J{T XX + 2T tx + T tt )(T xx - 2T tx + T tt 
Tee 



• 2 ' 

sin x 
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r/cr 



(Ttt - T xx + 3) 



(S + T XX + 4T^ 



T tt ) 



5-T, 



YX 



T 



I A 



m x 

[- T xx+ T xx(- S + ^- 2T ") 



XX ' XX 

+T XX (Tjf, (S + 2T„) -T« (H + Tu) - 4T 2 X ) 
+T V , (T« (S + T«) + 4T t 2 x ) -2T t 2 x S] / [65] , 

(101) 



By itself, (101) is not profound: the stress-tensor has 



4 independent components, and we have simply mapped 
them to a new set of variables (e, P,v,rja). The cru- 
cial question though is whether these variables behave as 
those of a conformal fluid. The answer appears to be yes, 
even for these initially highly distorted black holes 13 . 

To support this claim, we first note that in all cases 
we have l ooke d at to date, we are able to perform the 
inversion ( 101 1 to obtain real- valued hydrodynamic vari- 
ables that satisfy e > 0, P > and v E (—1,1). Fig. 
17 shows an example of v(t,\), from the r/, = 5 quasi- 
normal ringdown case discussed before. This suggests 
the stress tensor is consistent with that of some fluid. 
To check that it is a conformal fluid, in the top panel of 
Fig. 18 we show the quantity (e — 3P)/e resulting from 



the evolution at several different resolutions, using the 
same initial data (to contrast with Fig. Il3l here we are 



13 A few of the explicit examples shown here correspond to a 
"moderately" distorted initial black hole, with = 5 and 
c eq /c p (t = 0) » 1.2— sec Fig. [7] though the same conclusions 
hold for the most distorted cases we have looked at to date 
(c e g/cp(t = 0) ~ 2.2), w here the fluid velocity reaches a maxi- 
mum \v\ S3 0.54 — see Fig.[l7| 



JL iff 

2 4 

X 



n 




X 



FIG. 17: The velocity field v of the fluid describing the 
boundary CFT, displayed on a t, x spacetime diagram, ex- 
tracted from the = 5 quasi-normal black hole ringdown 
simulation described in Sec. |VI B| The first panel corre- 
sponds to the w y /w x —4 case, and the w y /w x = 32 case in 
the second panel is included for comparison. 



not merely testing that the stress tensor is traceless, but 
rather the more restrictive property that it arises from 
fluid flow where in the rest frame of the flow the pressure 
is isotropic) . At any given resolution, this quantity is cer- 
tainly nonzero due to truncation error; however, except 
for a transient at the initial time, the trends show that it 
is converging to zero, at or better than first order 14 . In 



That the rate of convergence is closer to first order is due to how 
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FIG. 18: The quantity (e — 3P)/e, which measures the ex- 
tent to which the hydro dynamic description is that of a 
fluid with equation of state e = 3P. The first panel is 
constructed from simulations run at 4 different resolutions, 
labeled by the mesh spacing relative to the highest resolu- 
tion run h/2. What is shown at each time is the L 2 -norm of 
(e— 3P)/e taken over the entire grid. Except possibly for an 
early transient (see Fig. 19 I, that this quantity is converg- 
ing to zero shows that e = 3P to within truncation error. 
In the second panel we interpret the finest resolution data 
as the continuum solution to within an uncertainty (shaded 
region), obtained from the magnitude of the difference be- 
tween the highest and next to highest resolution data (this 
assumes first order convergence of the extracted quantities). 



other words, this says that after the initial transient, the 
boundary stress tensor behaves like a fluid with equation 
of state e = 3P to within numerical truncation error. 
Fig. [i~9|is a close-up of the early transient behavior seen 
in Fig. |18[ and as before for the convergence of the so- 



lution discussed in Sec. V C it appears to be converging 
away in the sense that it affects a smaller region in time 
as resolution increases. 

An alternative way to present this information is dis- 



played in the bottom panels of Fig.'s |18| and [19| Here, 
assuming first order convergence of the extracted quan- 
tities, we take the finest resolution result as the contin- 
uum solution with an uncertainty (from truncation er- 
ror) given by the magnitude of the difference between 
the finest and next-to-finest resolution results; this error 
is shown as the shaded region about the finest resolution 
curve. For subsequent time-series plots we will display 



we extrapolate the solution to the boundary to extract the stress 
tensor components — the underlying solution for the metric still 
shows second-order convergence as discussed in Sec. IV CI 
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FIG. 19: An expanded view of Fig. [18] at early times. 



the data in this format. 

From this data we might also try to extract the lowest- 
order transport coefficient, namely the shear viscosity ?y, 
but here we encounter the two primary limitations of this 
sub-section's extraction method. First, this method only 
allows us to extract the product rja, so in order to find r\ 
we must divide by an independent calculation of a. Such 
a calculation may be achieved by substituting the ex- 
tracted velocity field into (96 1, which yields a |„, so we can 
compute r\a J a\ v w r\. Given our time-symmetric initial 
data, this o~\ v has zeros at t = 0, and periodically after 
that. Even though the "true" rja necessarily has zeros at 
exactly the same times, the extracted rja will not, because 
we have ignored the II^„ higher-order contributions in 
the map. Thus, such a calculation of rjo~/o~\ v ss r\ has the 
unattractive feature of diverging whenever o~\ v = 0. Sec- 
ondly, and more crucially, this method does not readily 
extend to allow the extraction of the higher-order trans- 
port coefficients. To remedy these shortcomings, we will 
now make use of a more general method to analyze the 
hydrodynamic behavior of the CFT T„ y . 



E. Higher-order Hydrodynamics 

In this section, we present an alternative method for 
comparing = (T Mi ,) CFT with the stress tensor of an 
J\f = 4 SYM fluid. To accomplish this, let us first add as 
many higher-order terms T\^ v to the hydrodynamic stress 
tensor (95) as is presently known, which includes viscous 



corrections up to second-order in the gradient expansion. 
This gives four additional higher-order transport coef- 
ficients to supplement the shear viscosity 77: the stress 
relaxation time tv, the shear vorticity coupling r u , the 
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FIG. 20: The mismatch AT^'* between the CFT stress 
tensor and the reconstructed hydrodynamic stress tensor 
with viscous corrections up to order i ma x = 0, 1, 2, for a 
simulation run with w y /w x — 4 initial data, whose final 
state black hole has horizon radius — 5. Here, ATgg lax = 

Toe — X^=cT Tgg , and the plots are normalized by the fluid 
energy density e. As discuss in the text and Fig. |18[ what is 
shown in each case is the data from the highest resolution 
run (dashed line) together with an estimated uncertainty 
(shaded region). Ideal hydrodynamics (top) shows a small 
residual mismatch; this is almost gone including first order 
viscous corrections (middle), and zero to within truncation 
error when second order corrections are added (bottom). 




FIG. 21: The same set of plots as described in Fig. [20] 
though here from the simulation of the w y /w x = 32 initial 
data. Notice the larger relative mismatches at a given order 



of the expansion compared to the w y /w x = 4 case in Fig. 20 
and that there is still a small residual even after including 
all corrections through second order. The trend going from 
ideal (top) to second order viscous hydrodynamics (bottom 
suggests that third and higher order corrections could fur 
ther reduce the mismatch. 



shear-shear coupling £ CT , and the Weyl curvature coupling 
£c- Our strategy here will be different from the one de- 
scribed in the previous sub-section, where we had first ex- 
tracted all hydrodynamic variables from the CFT stress 
tensor components, then exhibited evidence that the con- 
formal constitutive relations hold. Instead, we will now 
assume that the conformal constitutive relations hold at 
the outset, allowing us to reconstruct the hydrodynamic 
stress tensor order by order using only the energy density 
e and the fluid four- velocity u M . 



This reconstruction takes the form 

oo 

V = E T /$ (102) 

i=0 

where Tf$ corresponds to the stress tensor of an ideal rel- 

ativistic fluid, and Tj^J accounts for the i th -oxdet correc- 
tion in a velocity-gradient expansion. These corrections 
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are explicitly given in [55] for an M = 4 SYM fluid: 
T$ = eu IM u v + P± IMV 

1 



O'n P\p 



—a H (j a f3 



+ £,cC tla vf}U a U 



(103) 



where C^ av fi is the Weyl tensor and V is the Weyl co- 
variant derivative defined in [62] , One further input are 
the constitutive relations for our M — 4 SYM fluid: 

HV:2 (^T) 4 = 3P 



V = 



2 (-?r 



8tt 2 
Nc 2 
8tt 
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2ttT 
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2ttT 



(104) 



where N c 2 /(8tt) = l/(16-7rG), and T is the temperature 
of the fluid which we measure from the energy density e. 
The raw materials from which each T^t is built can be 
written in terms of CFT stress tensor components as 



f 



(T„ - T xx + E) 



-2Ti 



t-x 



-4T? X + (T U +T xx + Zf 



(105) 



where 5 is defined by ( 100 1 



We are now in the position to ask whether the bound- 
ary flow is consistent with an Af = 4 SYM fluid, and 
if so, the extent to which higher-order corrections are 
important in describing it. We address these questions 
by comparing each of the reconstructed stress tensors 
£W r C0 for imax = o, 1,2 to the full boundary CFT 
stress tensor of the numerical solution. The results of 
this comparison are summarized in Figs. [20] and [21] for 
a representative stress tensor component (normalized by 
the energy density e). The plots in Fig. 20 are from the 
Wy/wx = 4 solution, where the maximum velocity in the 
flow reaches ~ 0.f2 (see Fig. 17), and Fig. 21 is from 
the most asymmetric initial data case w y /w x — 32 that 
exhibits a maximum velocity of w 0.54. For the former 
case, the solution converges to behavior that is different 
from ideal hydrodynamics by at most sa 0.5% (modulo 
the transient near t = 0) ; the inclusion of first-order and 
second-order corrections successively decrease this differ- 
ence to a level which is essentially zero to within trunca- 
tion error. The analogous plots for the larger asymmetry 



case in Fig. 21 show deviations from ideal hydrodynamics 
of up to 2.5%; adding first order corrections reduces the 
maximum deviation to w 1.5%, and second-order correc- 
tions to just under « 1%. This is a rather remarkable 
level of consistency, in particular considering the trends 
in Fig. |2l] as succesively higher order viscous corrections 
are added. 



F. Passing to Minkowski Space 

Thus far, in the field theory dual of deformed black 
holes, we have focused on fluid flows on the boundary 
of global AdSs, namely I x S 3 . As is clear from sec- 
tion |VI C| and in particular from the second panel of 
Fig. |15| the fluid flow is essentially a compressive wave 
which starts with its peak at the equator of 5 3 and then 
travels to the poles and back, oscillating and damping 
out toward a static configuration where the temperature 
on 5 3 is everywhere constant. In this section, we would 
like to make closer contact with real-world fluid flows by 
conformally mapping our extracted boundary solution on 
K x 5 3 onto a corresponding flow in Minkowski space. 
To do this, we first explain the main concepts behind 
the conformal mapping; then we provide numerical re- 
sults based on our most anisotropic deformed black hole 
solutions. 

Up to a conformal transformation, M x S 3 is the cov- 
ering space of Minkowski space, R 3,1 . Therefore, we can 
obtain fluid flows on R 3 ' 1 by restricting our results to an 
appropriate patch of K x S 3 and then conformally map- 
ping to Minkowski space. There are many ways of doing 
this, simply because there are many ways of position- 
ing the Minkowski space patch within E x S 3 . We will 
mostly focus on one particular choice of patch which leads 
to a flow reminiscent of a head-on heavy ion collision — 
though, as on M. x S 3 , the initial conditions are charac- 
terized by full stopping rather than approximate rapidity 
independence. More specifically, the initial timeslice of 
our simulation, described as t = in K x S 3 , will cor- 
respond to a Minkowski timeslice t' = in which the 
fluid is stationary, and compressed into a region which is 
axisymmetric and significantly oblate. This is intended 
to be compared to the state of two heavy nuclei which 
have just achieved full overlap, though because the initial 
velocity profile is zero, a better analogy may be a quark- 
gluon plasma in a trap. The subsequent evolution of the 
fluid comprises two distinct types of expansion: longitu- 
dinal expansion along the axis of symmetry, and radial 
expansion. Overall, the flow preserves an 50(3) sub- 
group of the conformal group 5*0(4, 2) (as it must since 
the black hole to which it is dual has an 50(3) symme- 
try), but due to our choice of embedding of Minkowski 
space in K x 5 3 , this 50(3) is not the obvious one com- 
posed of rotations around a point in a single timeslice. 
Rather, it is the conformal 50(3) symmetry used in [63T 
156] to study generalizations of Bjorken flow. Rotations 
around the axis of symmetry of the Minkowski space flow 
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form an 5*0(2) subgroup of the conformal SO(3) symme- 
try. The rest of this SO (3) is composed of special confor- 
mal transformations, corresponding to conformal Killing 
vectors of Minkowski space. 

In order to map lx5 3 , covered by global coordinates 
x>i = (tiXi@i4>)i to IR 1,3 , covered by coordinates x a — 
(t',Xi,X2,X3), we use the transformations 
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X-2.IL 
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cost/L + cosx 

sin x 
(cost/L + cosx) 

sin x 
(cosi/L + cos x) 

sin x 
(cost/L + cosx) 



(106) 



The appearance of the AdS scale L on the right hand 
side of ( 106 1 is essential: only with this factor will ( 106 ) 
lead to a conformal mapping of the metric 



xS 3 
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= -dP 



dx^dx 
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sin 2 x 2 sin 2 8 2 



to the standard Minkowski metric 



d s 2 3il = fi b - dx a dx b = -(dt') 2 +(d Xl ) 



-(dx 2 ) 2 +(dx 3 ) 2 . 

(108) 



The appearance of L on the left hand side of ( 106 1 is 
inessential: it could be replaced by any quantity with 
dimensions of length. Doing so would amount to altering 
Minkowski space by an uniform dilation of both time and 



space. The conformal mapping ( 106 1 is accompanied by 

(109) 



the following rule for metric components: 

i dx^ 1 dx p Y c>3 
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dx a dx b ^ l 



9au 



where 



W 



cost/L + cosx 



(110) 



The Minkowski space patch on R x S 3 is the connected 
region including (£, x) = (0, 0) where W > 0. This region 
is easily seen to be the region — it < t/L < tt with < 
X <n-\t/L\. 

We have used coordinates (t,x,9,(f>) on I x S 3 , in- 
stead of our previous coordinates (t,x,8,4>), in order to 
preserve our freedom to position the patch in any way we 
wish on K x S 3 . For example, we could set t = t — to in 
order to "center" the Minkowski space patch on a global 
time t — t a . For our current purposes of describing a 
fluid flow reminiscent of a head-on heavy ion collision, 



the most useful choice is to set t = t, 



and 



y = arctan 
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sin x cos 9 



\J cos 2 x + sin 2 x sin 2 i 
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arctan (cot x esc ( 



(111) 



This mapping is an isometry of S 3 , and composing it 
with the mapping ( 106 ) gives 
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sin t/L 
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h sin x cos 9 




sinx 
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+ sin x cos 9) 




sinx 
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+ sin x cos 9) 




cosx 
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+ sin x cos 9) 



sin f cos i 



(112) 



Additionally, one may show by plugging (111 I into ( 110 ) 
that 



W 



1 



cos(f/L) + cos 9 sinx 



(113) 



The reason that the final mapping (112) is a good idea 
for producing a flow reminiscent of a central heavy ion 
collision is that the origin of Minkowski space maps to 
an equatorial point = (tt/2, 0) on S 3 , where the 

fluid density is at its peak on the initial timeslice. (Note 
that the coordinate <j> on S 3 becomes the usual angle <f> 
around the symmetry axis in Minkowski space, but the 
coordinate 9 on S 3 is more closely related to transverse 
radius from the symmetry axis in Minkowski space than 
it is to the angle of latitude with respect to the symmetry 
axis.) 

We now show how the boundary stress tensor Tj^* s = 
= (T Alzy ) CFT transforms under the above conformal 
mapping. First recall that we have effectively subtracted 
a term t M „ = hm^^'t^ from T^* s , by removing the 

term from the quasi- local stress tensor (21) prior 

to taking the boundary limit. Due to this subtraction, 
the correct rule for obtaining the stress tensor T* h ' on 
Minkowski space obtained via the conformal mapping de- 
scribed above is: 



T^ 1 = W- 



, dx» dx v K> 
dx a dx b » v 



(114) 



Had we not subtracted t^ v in our calculation of 



we would have had to include an additional term in ( 114 ), 
identified in [57] as a consequence of the holographic Weyl 
anomaly. In other words, the full stress tensor is not a co- 
variant object under the family of bulk diffcomorphisms 
that induce conformal transformations on the boundary. 
On the CFT, this anomalous term is not surprising: it is 
simply the usual Schwarzian derivative term that appears 
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FIG. 22: Temperature T denned in (115 1, for a simulation run with w y /w x = 32 initial 
data, whose final state black hole has horizon radius = 5. Each plot depicts the spatial 
profile of temperature in the Xi — X3 plane, taken at a constant t slice and at <j> = 0. One 
can recover the full spatial dependence by simply rotating each of these x\ — x$ profiles 
about the X3 axis. By interpreting the X3 axis as the beam-line direction, and xi as the 
radius in the transverse plane, the initial data in the first panel can be thought of as 
approximating a head-on heavy ion collision at its moment of impact. 



when one passes the stress tensor of an even-dimensional 
CFT through a conforms! transformation. From the bulk 
perspective, the anomaly originates from picking out a 
particular foliation before taking the boundary limit ( 19 1. 



To explicitly see that the rule (114) results in the 



correct stress tensor on Minkowski space, it is instruc- 
tive to consider the unsubtracted boundary stress ten- 
sor T®. = lim \( q 'T?,„. Since the anomalous term does 

not depend on the bulk dynamical variables g^, v , H^, </>, 
it suffices to evaluate this object in the case where 
§fiv = 0, Hfj, = 0,0 = 0. Performing this calculation in 
global AdSs, we use the foliator q = 1 — p = l/(l + r) and 
find that T^ xS — t^ u . On the other hand, this same 
calculation on the Poincare patch involves the foliator 
9 = + r 2 / L 2 ) cos(t / L) + r/Lsmxcos9) (essen- 

tially the Poincare radial coordinate) and one finds that 
= 0- Tnis calculation concretely demonstrates 
that T® u is not a covariant object. It additionally shows 
that by subtracting t^ from the boundary stress tensor 
on I x S 3 , one may pass the rest of the stress t ensor 



through the tensorial rule (114) 



and obtain the correct stress tensor on 



to have an in- 



If we now take the stress tensor T^* s 

viscid hydrodynamical form T^* s — eu^u^ + e /3 
as was demonstrated in Sec. |VI E[ then one can straight - 

^ and e 83 ' 1 = 



forwardly show that " 

T/ j/-4 £ Rxs 3 Tjp to a cons tant factor specifying the num- 
ber of degrees of freedom, the temperature of a conformal 
fluid is T = e 1 / 4 . Setting L — 1, the quantity we are go- 
ing to plot is 



where T = T R ' is (up to the aforementioned con- 
stant factor) the temperature in Minkowski space and 



lxS J 



is the energy density on K x S 3 . 



W 



-l e l/4 



We are justified in plotting temperature starting at 
Minkowski time t' = because, by construction, our ini- 
tial data perfectly conforms to the inviscid hydrodynamic 
ansatz. Moreover, again by construction, the velocity 
field vanishes exactly on the initial timeslice. Thus we 
are describing a version of the Landau-Khalatnikov full- 
stopping scenario [68j [69] , but with radial flow explic- 
itly included. The results of section |VI E| in particular 
Figs. [20] and [21] show that hydrodynamics remains ap- 
proximately valid throughout the evolution, with modest 
contributions from shear viscosity and small contribu- 
tions from second-order transport coefficients. 

The four panels of Fig. [22] shows the temperature T 
in Minkowski space reconstructed from a simulation 
with initial asymmetry w y /w x = 32. Each panel in 
Fig. 22 is taken from at' — const slice, and at <f> = 
(i.e. X2 — 0). Since the transformation to Minkowski 
space preserves the original ^-symmetry in global 
AdS, one can recover the full spatial profile by simply 
rotating each of these x\ — X3 profiles about the X3 axis. 
To quantify the initial anisotropy in the temperature 
profile, we use ((^i) 2 ) / ((X3) 2 ), where (A(xi, 0, £3)) 
j dxidx3T(xi, 0, X3) i A(xi, 0, £3)/ / dxidx^T^xi, 0, X3) 4 
is the energy-weighted average of some function 
A(xi, 0,2:3) defined on the initial t' = slice and at 
= 0. This quantity evaluates to 6.94 for the data 
displayed in Fig. [22] We also calculate the widths 
<5i and 6 3 given by T 4 (£ l5 0,0) = T 4 (0,0,0)/2 and 
(115) T 4 (0,0,5 3 ) = T 4 (0,0,0)/2, then construct the ratio 
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(5i) 2 / (S^) 2 this evaluates to 18.7 with the data 



VII. CONCLUSION 

We have described a generalized harmonic scheme for 
solving the Einstein field equations on asymptotically 
anti-de Sitter spacetimes in 4+1 dimensions. Though 
restricting to SO(3) symmetry in this initial code, we ex- 
pect that many of the methods developed here to achieve 
stable, consistent evolution will carry over to scenarios 
with less symmetry. This will be needed to tackle the 
main motivation for developing this code, namely to ob- 
tain the gravity dual to a heavy ion collision, and its 
subsequent quark-gluon plasma formation. 

As a first application, we studied the quasi-normal 
ringdown of highly distorted black holes, and the cor- 
responding behavior of the stress energy tensor of the 
dual CFT on the boundary of the spacetime. We find 
quasi-normal mode frequencies that are consistent with 
previously published linear modes, as well as modes that 
can be modeled as arising from non-linear mode-coupling. 
We further find purely decaying modes that we attribute 
to gauge. To be certain this is the case would require 
transforming the solution to coordinates fully compati- 
ble with perturbative calculations, which is a rather non- 
trivial task numerically, and so we relegate this to a fu- 
ture study. 

The boundary stress energy tensor exhibits corre- 
spondingly large initial fluctuations, yet we find that its 
behavior is consistent to better than 1% with that of 
an TV = 4 SYM fluid with equation of state e = 3P, 
and with corresponding transport coefficients, essentially 
from t = onwards. This is in contrast to the numeri- 
cal results reported in ! 23-26] , where only after a certain 
time did the boundary behavior approach that of a fluid. 
However, those studies looked at scenarios more akin to 
black hole formation, i.e. beginning with states that do 
not contain large black holes (or black branes), with the 
black holes forming at later times. Such processes are ex- 
pected to be dual to thermalization, whereas our study is 
that of equilibration beginning from an inhomogeneous 
though thermal state. Nevertheless, it is curious that 
the link between Einstein and Navier-Stokes seems to be 
holding even in these far-from-equilibrium scenarios, and 
it will be interesting in future work to explore how far 
this relationship extends. Recently it was shown that 
the Rayleigh-Plateau instability in a fluid stream and 
the Gregory-Laflamme instability of a black string are at 
least qualitiatively similar |70j . These findings similarly 
suggest that, in some situations, the physics described by 



the Einstein and Navier-Stokes equations could exhibit 
similarities even in the most non-linear, near-singular 
regimes. 

By passing to an appropriate Minkowski patch of the 
boundary of global AdSs , we are able to extract a fluid 
flow which starts from a compressed disk and leads to 
expansion both in the radial and longitudinal directions. 
The SO(3) symmetry is still present, but as a conformal 
symmetry rather than an isometry of the flow. This flow 
exhibits the main qualitative features of the full-stopping 
scenario in heavy ion collisions, which, though gener- 
ally considered implausible in light of asymptotic free- 
dom, has some interesting phenomenological successes 
[71] . It would be interesting to pass our flow through 
a hadronization algorithm and extract the rapidity dis- 
tribution of produced particles. We leave this task for 
future work. 

A further direction for future study, that could be ac- 
complished within the simplification of the current SO '(3) 
symmetry, is to couple gravity to additional forms of 
matter. In the context of AdS/CFT, it is interesting 
to consider tachyonic scalar fields, since these amount 
to relevant operator insertions in the boundary CFT. As 
one might expect, the resulting deformation is dramatic: 
on the gravity side, these tachyonic scalars back-react in 
such a way as to significantly change the boundary con- 
ditions of the metric. The evolution of such systems in 
the full non-linear regime is unknown, though hopefully 
methods similar to those introduced here will be able to 
handle such spacetimes with "deformed" AAdS bound- 
aries. These methods are also applicable to AAdS space- 
times with different boundary topology, which (with the 
appropriate matter content) are of interest to condensed 
matter applications of the duality. 
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The extent to which (Oi) 2 ) / (O3) 2 ) and (<5l ) 2 /{<5 3 ) 2 differ re- 
flects the extent to which the temperature profile deviates from a 
perfect Gaussian. As a measure of anisotropy, the ratio of widths 
at half-max (<5i) 2 /(S3) 2 is the more faithful of the two. 
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Appendix A: Effect of Scalar Backreaction on 
Metric Fall-off 



Here we compute the effect that a static spherical dis- 
tribution of massless scalar would have on the AAdSs 
metric. We follow the discussion in [35]. The configura- 
tion only has radial dependence <j) = <f>(r), so the metric 
must take the form 



ds 2 = - (l + r 2 + 0{r- x j) dt 2 + 



1 



dr 2 



(Al) 



where (x(r) grows slower than r 3 as r — > oo in order to 
preserve the value of the cosmological constant. We have 
set L = 1, so that A 5 = —6. In a previous section, we 
had used the AdSs Klein-Gordon equation ( 16 ) in spher- 



ical symmetry to find the general leading order behav- 
ior (17) of scalar fields in AdSs. We now find the metric 



fall-off behavior when scalar back-reaction is taken into 
account, limiting ourselves to the h pp component of the 
metric deviation and using the compactified coordinate 
p = r/(l+r). In order to find this fall-off, the strategy 
is to solve for p,(r) as a power series in r, from which we 
can reconstruct 



hpp — 



MM 



(A2) 



where we have used dr/dp = 1/(1 — p) 2 ~ r 2 along with 



g rr and 



(1 



p{r)/r) 



p(r) 



+ 0(r~ 5 ) 



(A3) 



We can deduce the asymptotic behavior of p(r) from 
the the Hamiltonian constraint 16 



(4) i?-2A 5 = 1%tt Pe . 



(A4) 



where 



(r + r 3 ) 2 



(r + r 3 ) 2 



(A5) 

is the Ricci scalar associated with the 4-metric on the 
slice, and 

PE = \(l + r 2 - M (r)/r)(0'(r)) 2 + \m 2 <? (A6) 

is the scalar field energy density. Keeping only the terms 
in (A4| that dominate at large r, we obtain the form of 



the Hamiltonian constraint near the boundary 
p'(r) p(r) 



r 2 ((t>'{r)) 2 +m 2 <i> 2 . (A7) 

The scalar field goes as (p(r) ~ r~ A near the boundary 
for some A, so we see that the right-hand-side of (A7) 
goes as r~ 2A . Matching the left and right-hand sides, we 
conclude that 



p(r) 



and consequently, 



„3-2A 



„2-2A 



(A8) 



(A9) 



To make use of what we have learnt, first sup pose 
that we have a scalar with negative m 2 < 0. From (17) 



and (18), a non-zero A branch describes the leading be- 
havior for large r, and is defined by a fall- off pow er of 
A = A + = 2 - V4 + m 2 . Then equation ( A9) pre- 



dicts that matter backreaction induces a metric devia- 
tion hpp ~ r -2+2 v / 4+m 2 fagfc jg larger than the vacuum 



metric fall-off we assumed in (54) and in (55). This case 
corresponds to a tachyonic scalar configuration. These 
configurations are stable in AdS as long as its mass sat- 
isfies the Breitenlohner-Freedman bound |72] given by 
-(£) - 1) 2 /(4L 2 ) = -4 < to 2 in D = 5 dimensions: it is 
stable in the sense that its conserved energy functional 
remains positive in the function space of all fluctuations 
for which the energy functional remains a convergent in- 
tegral. The qualitative reason for the positivity of this 
energy functional in AdS is that despite the negative un- 
bounded potential of the tachyonic scalar, the positive 
kinetic terms in the energy functional dominate the neg- 
ative potential terms so long as the scalar falls off suffi- 
ciently quickly near the boundary 17 . We shall elaborate 
on the behavior of these fields in a subsequent paper. 
For the present paper, we consider the case of scalars 



with m > 0. Again looking at (17), (18), notice that 



we must now turn off the A branch in order to have a 
scalar field that vanishes at infinity, which is required 
for any scalar field configuration with finite energy. The 
result is a localized matter distribution defined by a fall- 
off power A = A + = 2 + \/4 + ra 2 . Applying our result 
(A9) to this case, we find that backreaction induces a 
metric deviation h pp ~ r -2-2\/4+m 2 fas,t is subleading 
compared to the vacuum metric fall-off (55), and thus 
leaving it unchanged. 



A scalar with the critical fall-off behavior goes as 



for 



See Sec. 



Ill A 



A c = (D - l)/2 = 2 in D = 5 dimensions; 
inferred from this observation. 



the BF bound can be 
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Appendix B: Boundary Conditions for Linearized 
Gravitational Perturbations of AdS 



Fefferman- Graham Theorem 

A theorem due to Fefferman and Graham [36] states 
that a distinguished coordinate system (z, x l ) exists near 
the boundary of any asymptotically AdSs spacetime M 
in which the metric takes the form 



L 



G = „ [dz 2 + gijdx l dx j ) 



(Bl) 



and that there is a convergent power series solution for 
the gij coefficients in (Bl) given by 



9ij = 9(o)ij + 9(2)ijZ 2 + 9U)ijZ + 2h(4,)ijZ log z + 0(z 6 ) 

(B2) 

(the boundary dM is located at z — 0). 

Let us first translate this statement in terms of global 
AdS 5 coordinates (t, r, x, Q)4>)- We follow the discussion 
in [67j . though we use a different Fefferman- Graham ra- 
dial coordinate z from the one use there, whose relation 
with the global radial coordinate r is given by 



1 



(B3) 



To leading order, this relation between z and r near the 
boundary at z = 0, or equivalently r — > oo, simply reads 



1 



0(l/r 2 ). 



(B4) 



In terms of the global coordinate q = 1/(1 + r) in- 
troduced in Sec. II E to leading-order we have z — 
1/r + 0(l/r 2 ) = q + 0(q 2 ) near the q = boundary. 
So the Fefferman-Graham theorem's power series solu- 



tion in these coordinates reads just like (B2 ), replacing z 
by q 



9ij = 9{0)ij + 9(2)ij<f + 9(i)i 3 <l A + 1h 



^m^Ogq + Oiq 6 ) 
(B5) 

0). 



(the boundary dM is located at q 

The 9(o)ij corresponds to the non-radial components 
of the pure AdS$ metric. The g(2)ij can be expressed, 
as was done explicitly in [37] , in terms of the Ricci ten- 
sor R(o)ij and the scalar curvature i?( ) constructed from 
9(o)ij, and is thus fixed. The leading-order dynamics then 
first appears 18 in the terms controlled by g(A)ij- Inspect- 
ing (B5|, this implies that the leading-order asymptotics 
is Gij ~ q 2 , in agreement with the boundary conditions 
we wrote down in Sec. IIVBI 



Ishibashi-Wald Boundary Conditions 

Here we show how our metric boundary conditions re- 
late to those derived by Ishibashi and Wald in [34] . Work- 
ing perturbatively, they derive conditions under which a 
metric perturbation must have its boundary behavior en- 
forced by explicit boundary conditions. The fields gov- 
erning the perturbations, which we collectively denote 
by /, are decomposed in terms of spherical harmonics 
§£;(£!*) according to 



1 



(B6) 



The $ fields for the tensor and vectors modes of the grav- 
itational perturbations are shown to require no boundary 
conditions at infinity (see Proposition 3.1 and its preced- 
ing discussion in [34]). However, the scalar modes of the 
gravitational perturbations require boundary conditions, 
and their corresponding $ are shown to behave asymp- 
totically as 



$ = F {r) [2a log(l/r) + b + 2L 2 /r 2 log(l/r) + 
where the function Fo(r) is asymptotically given by 



(B7) 



Fa(r) 



1 

h72' 



(B8) 



In terms of the global coordinate q, the implication is 
that for the scalar modes of the gravitational perturba- 
tions the near-boundary behavior is 



G, 



2a q 2 log(g) + b Q q 2 + 2c q 4 log(g) 



(B9) 



For the scalar modes, imposing boundary conditions 
amounts to setting the ratio bo/ao. Comparing (B9) 



to (Bl ),(B5), we see that the choice picked out by coordi- 



nates that are Fefferman-Graham-like near the boundary 
is bo/ao — > oo i.e. ao = 0. This removes the logarithmic 
branch controlled by ao, so the leading-order asymptotics 
is Gij ~ q 2 , as desired. 



In D = d + 1 dimensions, the leading-order dynamics appears in 
the 9(tj)ij term, since gr n )ij for n < d are determined by 9(o)ij- 
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Appendix C: Tables of Scalar Quasinormal Mode 
Frequencies 



fund. 


k=0 








k= 


=2 








k=4 




Th 

L 


=12.2 


(2.95 


± 


0.13)2 - 


«(2.67±0.02)^ 


(3 


96 


± 


0.03)^ - 


i(2.33dz0.07)^ 


(4.93 zb 0.007)^ 


-«(1.83± 0.05)2 


r h 
L 


=11.3 


(2.95 


± 


0.14)2 - 


i(2.64dz0.04)^ 


(3 


94 


± 


0.02)^ - 


i(2.30± 0.08)2 


(4.90 zb 0.003)2 


-i(1.80zb 0.04)2 


rh 
L 


=10.5 


(2.96 


± 


0.14)2 - 


i(2.64zb0.05)^ 


(3 


95 


± 


0.01)2- 


i(2.29±0.08)^ 


(4.92 zb 0.002)^ 


-i(1.78± 0.04)2 


r h 
L 


=9.0 


(2.99 


zb 


0.15)2 - 


i(2.63zb 0.07)2 


(3 
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dz 
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Th 

L 


=6.5 
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zb 


0.14)^ - 
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(3 


96 


dz 
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(4.80 zb 0.01)2 - 


-i(1.65zb 0.003)^ 


rn 
L 


=5.0 


(3.11 


dz 


0.09)^ - 


i(2.62±0.20)^ 


(3 


94 


dz 


0.05)^ - 


i(2.21dz0.04)^ 


(4.79 ±0.08)^ - 


-i(1.60zb0.02)a 


r h 
L 


=4.5 


(3.18 


± 


0.06)2 - 


i(2.67±0.22)^ 


(4 


00 


± 
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i(2.24±0.04)^ 
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-i(1.69zb 0.04)2 


Th 

L 


=4.0 


(3.21 


± 
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-i(2.67±0.23)^ 


(3 


99 


dz 
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-i(1.76zb 0.04)2 


r h 
T, 


=3.3 


(3.33 


dz 


0.08)^ - 


i(2.71zb0.20)2 


(4 


05 


± 


0.08)^ - 


i(2.27± 0.03)^ 


(4.9 zb 0.1)2 -i 


(1.89 zb 0.05)2 



TABLE IV: The fundamental (n = 0) quasi-normal mode frequencies ui r — ioji ext racted 
from the scalar field variable cb, from evolution of initial data as described in Sec. VI B 



These are shown for various horizon radii ru of the final state AdS-Schwarzschild black 
hole with 50(4) quantum numbers k > 0, 1 = 0, m = 0. Uncertainties are estimated from 
convergence studies. 
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TABLE V: The first overtone (n = 1) quasi-normal mode frequencies ui r — iuui ext racted 
from the scalar field variable 6, from evolution of initial data as described in Sec. VI B 



These are shown for various horizon radii of the final state AdS-Schwarzschild black 
hole and for harmonics with 5*0(4) quantum numbers k > 0, / = 0, m = 0. Uncertainties 
are estimated from convergence studies. 
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